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Polynomial  Approximation  Techniques  for  Differential  Equations 

in  Electrochemical  Problems 


B.  Stanley  Pons 
Department  of  Chemistry 
University  of  Alberta 
Edmonton,  Alberta,  Canada 
T6G  2G2 


I .  INTRODUCTION 


Stephen  Feldberg's  finite  difference  methods  [1]  for  digital 
simulation  has  been  the  primary  method  of  approximating  the  solutions 
to  second  order  partial  differential  equations  found  in  electro¬ 
chemistry.  The  intensity  of  usage  of  digital  simulations  by  these 
finite  difference  methods  has  increased  markedly,  and  it  has  become 
apparent  in  some  cases  that  there  is  a  need  for  faster  techniques 
capable  of  delivering  accurate  results  as  well  as  being  able  to 
simulate  extremely  "fast"  systems  more  efficiently.  As  a  result, 
several  groups  [2]  have  made  significant  contributions  toward  im¬ 
proving  finite  difference  methods  for  electrochemical  problems. 

But  instability  and  lengthy  computational  times  under  certain  con¬ 
ditions  still  manifests  itself  as  a  primary  drawback. 

The  methods  described  in  this  chapter  provide  some  of  the 
answers  to  these  problems.  The  mathematics  involved  in  these 
polynomial  approximation  methods  are  only  slightly  more  com¬ 
plicated  than  the  finite  difference  equations,  but  the  accuracy  of 
the  results,  stability  limits,  and  speed  of  execution  of  the  com¬ 
puter  programs  make  their  use  well  worth  the  small  extra  effort 
needed  to  implement  them. 

It  is  not  intended  to  imply  that  the  methods  described  are 
the  "ultimate"  electrochemical  differential  equation  simulations. 

The  rapid  convergence  of  low  approximation  order  solutions  to  the 
exact  solutions  has  been  proven  mathematically  for  many  types  of 
problems,  but  certainly  not  all.  However,  it  is  this  type  of  con¬ 
vergence  in  a  method  that  urges  us  to  seek  general  theorems  that 
may  lead  to  faster,  even  more  effective  methods  of  approximation. 
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The  interest  in  the  use  of  new  highly  efficient  numerical 
methods  for  the  solution  of  complex  mathematical  problems  is  further 
evidenced  by  the  overwhelming  number  of  publications  on  the  sub¬ 
ject  that  have  appeared  in  the  recent  mathematics  and  engineering 
literature  [3] .  The  success  of  these  techniques  may  be  attributed 
to:  a)  the  impressive  degree  of  accuracy  that  is  obtainable; 

b)  the  fact  that  many  of  the  methods  may  be  constructed 
from  only  a  few  basic  principles; 

c)  the  consideration  that  many  methods  are  "modular"  in 
construction,  and  the  "modules"  are  available  as  highly 
efficient  algorithms  amenable  to  digital  computer 
implementation . 

Harrison  and  Gray  [4]  used  Chebyshev  polynomial  approximations 
some  years  ago  for  an  electrochemical  simulation.  More  recently, 
Whiting  and  Carr  [5]  introduced  an  advanced  and  highly  efficient 
method  (orthogonal  collocation)  into  the  simulation  of  electrode 
reactions.  Their  work  outlines  the  theory  of  orthogonal  colloc¬ 
ation  and  its  application  to  several  electrochemical  mechanisms 
occurring  during  chronoamperometric  experiments.  Bewick,  Mellor, 
and  Pons  [6]  applied  the  technique  to  some  actual  complicated 
systems  and  extended  the  method  to  simulate  modulated  specular 
reflectance  transient  responses  of  the  intermediates  formed  during 
the  reaction.  More  recently,  Rieker  and  Speiser  [7]  have  used 
orthogonal  collocation  to  simulate  a  full  spectrum  of  cyclic  volt- 
ammetric  responses. 

The  methods  that  will  be  discussed  in  this  chapter  are 
suitable  in  nature  to  the  solution  of  a  great  many  types  of  dif¬ 
ferential  equations.  As  Whiting  and  Carr  (5)  pointed  out,  it  is 
usually  just  a  matter  of  changing  a  few  steps  in  the  program  to 
switch  from  one  problem  to  one  of  a  completely  different  nature, 
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e.g.,  the  chronoamperomecric  response  of  an  e.c.e.  mechanism  in 
electrochemistry  to  a  dynamic  model  of  a  differential  scanning 
calorimeter.  With  this  in  mind,  it  is  also  important  to  remember 
that  even  highly  efficient  methods  may  be  improved  by  using  cer¬ 
tain  techniques  to  treat  peculiar  situations  that  arise  in  specific 
problems.  I  have,  however,  tried  to  keep  the  methods  as  general 
as  possible. 

Special  attention  has  been  given  to  the  integration  of 
"stiff"  differential  equations,  since  it  is  with  this  type  that 
electrochemists  are  most  often  confronted  in  their  routine  work. 
Equations  of  this  type  are  of  considerable  interest  in  engineering 
research  presently,  and  new,  more  efficient  methods  are  being 
developed  rapidly. 

Some  simple  approximation  methods  are  also  described  that 
may  be  solved  on  a  calculator.  These  methods  are  of  acceptable 
accuracy  in  many  types  of  experiments,  and  some  are  of  accuracy 
approaching  the  more  sophisticated  methods. 

II.  THE  METHOD  OF  WEIGHTED  RESIDUALS 
A.  Introduction 

Consider  the  boundary  value  problem 

0t=0xx  <  1) 


where 


a  -  60(x,t) 

*t  ~  fit 


(  2) 


4 


and 


0 


xx 


6x2 


3) 


We  assume  that  0  is  defined  and  continuous  in  the  domain  W.  We 
also  define  0q(x)  -  0(x,O)  as  the  initial  condition  of  0,  and 
0z(x,t)  are  the  values  of  0  at  the  boundaries  Z  of  the  domain. 

We  choose  an  approximate  solution,  0 A  to  the  problem  in  the  form: 


n 


0h  =  0z=o(x,t)+0z=L(x,t)+£  ^oi(t)0i(x) 


4) 


where  0A  is  the  approximate  solution,  0z=q  an(*  0Z_L  are  t*ie 
solutions  at  the  boundaries  0  and  L,  and  the  0 ^  are  the  basis 
functions  (which  may  be  polynomial,  trigonometric,  or  other  types). 
Various  techniques  would  prescribe  whether  the  basis  functions 
satisfy  the  differential  equation,  the  boundary  conditions,  or 
both. 

The  residual  of  equation  (1)  is  defined  as 


R<0)  -  0XX  -  0t  (5) 

so  that  if  0A  is  an  accurate  approximation  to  the  actual  0,  then 
R(0a)  will  equal  zero.  It  is  the  purpose  in  the  weighted  residual 
methods  to  pick  the  (t)  such  that  the  residual  tends  to  zero. 

The  a^(t)  are  chosen  by  specifying  that  the  integral  of  the  re¬ 
sidual  times  some  weighting  function  w  is  equal  to  zero.  Thus  it 
is  the  "weighted  average"  of  the  residual  that  is  specified  to  be 


zero  over  the  whole  domain: 


/w.R(fL)dx  =0  (6) 

W  3  A 

The  choice  of  the  weighting  function  w^  determines  which  method 
of  weighted  residuals  is  being  utilized.  Some  of  the  more  pop¬ 
ular  methods  are  listed  here  [8]  . 

1.  Integral  Method  (method  of  moments,  1st  order  subdomain) 

* 

The  Wj  are  chosen  to  be  x-3,  j=0,l,2,3, -  So  for  the  first 

order  approximation  problem,  x-^l,  the  integral 

/R(0.)dx  =0  (7) 

W  A 

is  evaluated.  The  result,  for  time  dependent  cu,  is  a  first 
order  differential  equation  in  0 A. 

2 .  Galerkin  Method 

In  this,  a  most  widely  used  and  highly  accurate  method, 
the  basis  functions  for  the  0 A  are  chosen  as  the  weighting 
function,  i.e., 

w.  =  0.  (x)  (  8) 

J  J 

Because  of  this  choice,  and  the  completeness  of  the  basis  set, 
the  method  can  be  made  exact  as  j  00 . 

3.  Variational  Methods 

Using  the  calculus  of  variations,  certain  problems  may  be 
approximated  similar  to  weighted  residual  methods.  The 
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solutions,  when  accessible,  compare  in  accuracy  to  the  Galerkin 
method.  However,  the  method  is  not  directly  applicable  to  all 
general  second  order  partial  differential  equations,  and  is 
mentioned  here  for  reference  only. 


4.  Least  Squares  Method 

Here,  the  weighting  functions  are  chosen  to  be 


<5R  <0A) 


j  ”  SoTTtT 


(  9) 


The  cu(t)  are  provided  by  n  equations  from  the  set  of  equations 


An 

SoTTtT  ■  2rR(0^hT7m  dx  ■  1  ■ 


(  10) 


for  a  one  dimensional  situation. 


5.  General  Collocation 


The  weighting  functions  are  given  by 


Wj  =  6  U-Xj)  , 


(  ID 


6  being  the  Dirac  delta  function,  which  is  defined  as  follows: 


6 (r)  =  0  ,  r^O 

/f(r)«(r)dT  =  f<0) 


(  12) 


Note  that  if  f(r)  =  1,  /6(r)dx  =  1.  Also,  it  is  true  that  6 (r) 
is  only  defined  during  an  integration. 
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It  is  then  immediately  obvious  that  this  choice  of  we¬ 
ighting  function  leads  to 


0 


/w.R(jZL)dx  =  /6  (x-x,)R(0.)dx  =  R(0,(x.)) 
W  3  A  D  A  A  J 


(  13) 


so  that  the  residual  is  given  in  terms  of  a  fixed  set  of  x_ 


only.  The  x^  are  any  set  of  points  in  the  integration  interval. 


6 .  Orthogonal  Collocation  Method 

This  highly  accurate  method,  which  is  used  in  this  work 
for  accurate  simulations  of  complicated  systems,  is  the  same 
as  method  (5)  above  with  the  collocation  points  x^  chosen  as 
the  real  roots  of  an  orthogonal  polynomial. 


7 .  Subdomain  Method 

The  integration  domain  W  is  subdivided  into  subdomains  W/n, 

n=l,2,3, - .  The  weighting  functions  are  equal  to  1  when  x  is 

in  a  particular  subdomain,  and  equal  to  zero  when  it  is  not. 

If  the  number  of  subdomains  is  increased,  the  residual  ap¬ 
proximates  zero  closely  in  more  regions,  and  eventually  ap¬ 
proximates  zero  over  the  entire  region. 


In  the  following  section,  a  very  useful  approximation  is  made 
to  the  general  diffusion  problem.  Several  weighted  residual  me¬ 
thods  are  then  used  to  provide  the  concentration  profiles,  and  a 
comparison  is  made  of  the  results. 


B.  The  Diffusion  Boundary  Layer  Approximation 

As  an  example  of  one  method  that  may  be  used  on  a  hand  calc¬ 
ulator,  we  develop  an  approximation  for  the  EC  mechanism  using  a 
boundary  layer  approach.  We  have: 
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at  a  planar  electrode  under  conditions  of  semi  infinite  linear  di 
f fusion,  and  the  conditions: 


[A] 


X,0 


[A°] 


(  15) 


[B] 


X,0 


[A] 0/T 


0 


(  16) 


(  17) 


This  problem  is  one  of  a  potential  step  applied  to  the  elec¬ 
trode  of  such  magnitude  that  the  surface  concentration  of  species 
A  is  reduced  immediately  to  zero.  We  wish  to  determine  the  con¬ 
centration  profile  of  species  B  as  a  function  of  time.  The  dif¬ 
ferential  equation  to  be  solved  is 


6[B] 

6T 


D  iilfL  -  mb] 

B  6X2 


(  18) 


Assume  that  L  is  some  distance  from  the  electrode  such  that  at  the 
end  of  the  experiment,  there  has  been  no  diffusion  of  any  species 
to  L.  In  this  case,  we  define  the  dimensionless  parameters 


t 


X  -  XL2 
L  B" 


C 


-  11 L 

B  [  A°  1 


(  19) 


Inserting  these  parameters  into  (18),  we  obtain 


5cb  {2cb 

6t  '  «x2 


(  20) 
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The  initial  and  boundary  conditions  are  also  changed  to 


cB(x,0)  =  cB(l,t)  =0  (21) 

cA(0,t)  =1  (22) 


Now  consider  the  diffusion  boundary  layer.  When  the  potential 
step  is  applied,  B  immediately  starts  diffusing  into  the  solution 
toward  the  cell  wall  which  we  assume  to  be  no  closer  than  L.  As¬ 
sume  that  the  concentration  of  B  ahead  of  this  boundary,  whose 
distance  from  the  electrode  we  shall  call  b,  is  equal  to  zero.  We 
also  assume  that  behind  the  boundary,  B  is  represented  by  a  con¬ 
tinuous  gradient.  We  note  that  b  is  a  function  of  t  only.  We 
may,  under  this  supposition,  treat  distance  from  the  electrode 
surface  as  some  fraction  of  the  parameter  b(t).  We  now  define  a 
dimensionless  distance  transformation,  which  is  simply  a  fraction 
of  the  diffusion  boundary  thickness: 


Since  the  concentration  is  to  be  expressed  as  a  function  of  distance 
from  the  electrode  surface,  we  have 

CB  =  cB(u)  =  ^(bltf)  (  24) 

Note  that  we  have  combined  distance  and  time  into  one  new  in¬ 
dependent  variable  p. 
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From  the  nature  of  the  model,  the  following  general  boundary 

conditions  must  be  imposed: 

(1)  The  dimensionless  concentration  at  the  electrode  surface 
is  1,  i.e. 

cB  =  f(y)=f(0)=l  (25) 

(2)  At  the  edge  of  the  moving  boundary  (y=l) ,  the  concentra¬ 
tion  is  0,  i.e. 

cB  =  f(y)=f(l)=0  (26) 

(3)  There  is  no  flux  of  material  across  the  moving  boundary,  i.e. 

f’(l)=0  (27) 

(4)  The  moving  boundary  initially  is  located  at  the  electrode 
surface,  i.e. 


b (0) =0 


(28) 


Now  we  rewrite  the  differential  equation  in  terms  of  the  new 
single  independent  variable  y. 


cb  =  f^>“f*bTtT) 


(29) 


6c 


5 1 


B  6  f  (y  ) 
6 1 


x 


X 


-) 


6 1 


6  ( 


=  _f  '  (_* - )  {_* - )  ( 

'  U  /  i-  \  •  '  /  i_  \  *  \ 


b(t) 
1 


) 


b  (t)  'b(t)'  lb(t) 


(— — ) 
'  U  t  *-  \  * 


...  6  f  b  ( t)  '  _  if  bit)  7  6  b  ( t) 


) 


fit 


db  ( t) 
dt 


=  =f ’  (y)  y  b(t)-1b’ (t) 


(30) 
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6cB_6f(FTtr)  _ef<F(tT)  6(b(t))  -MljW  1 
sir  &  . Tx  :  6x  f  (p)(b(tr>,B£  (w)b  (t) 

6  (BTtr) 


y- 6  (t,<^’b~1(t)  ^TiT^-g Wt,,fc)  'eW1  •  f  (v)b(t)-2  (32) 
{x  6(FTtT 


The  differential  equation  (20)  now  reads: 


-f  *  <VJ)  V  MtrVtt)  -  b“2f"(y)  -  &f(y) 


The  residual  of  (33)  is  simply  the  left  hand  side  minus  the  right 
hand  side.  Now  we  choose  an  approximating  function  for  f(y).  The 
simplist  polynomial  approximation  that  satisfies  all  of  the  boundary 
conditions  (25-28)  is: 


f(y)  -  (1-y) 2  «  l-2u  +  v2 
We  have: 

f* (w)  «  2y  -  2 


f"  (y) 


Rj.  *  -y(2y-2)b'(t)b(t)’1  -  2b(t)‘2  +  8(l-2p+y2) 


or  b(t)2R£-  (2vi-2y2)b'(t)b(t)  -  2  +  &b2lt)  (l-2v+p2) 


Now  we  shall  apply  various  MWR  to  solve  for  the  moving  boundary 
position  b.  Insertion  of  this  value  into  equation  (34)  will  give 


the  concentration  profiles  for  species  B.  The  results  for  each 
method  are  then  compared  to  the  exact  solution. 
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Method  (1) :  The  integral  method  (w^*l) 

For  the  highest  accuracy,  we  demand  that  the  residual  be  equal 
to  sero  in  the  weighted  average  sense: 

bVw.Rgdy  «  {lR£dy  *  0 

i.e. 

!'  (2y-2p2)b’  (t)b(t)  -  2  +  gb (t) 2  (l-2p+vj2)dy  *  0 
0 

/>  -  feMtiblt i  .  2  +  .  0 

The  first  order  differential  equation  is  solved  for  b(t): 

b(t)b*  (t)  -  6  +  6b2 (t)  -  0  (42) 

and  since  b'  (t)  «  ,  we  have 


(39) 

(40) 

(41) 


b(t)db(t) 

(8b2-6) 


dt 


(43) 


which  is  integrated  to  give 

Tc  1 1/2 

b(t)  -  |(l-exp(-2gt)>  (44) 

This  approximation  describes  the  movement  of  the  diffusion  boundary 
layer  in  time. 

The  approximate  concentration  change  with  distance  and  time  is 
obtained  by  substitution  of  (44)  into  (34): 
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cB  -  f <v>  -  (1-v)2  -  (1-^y)2 


(X  '  f< 


exp(-2gt)) 


The  solution  is  valid  until  the  moving  diffusion  boundary  layer 
reaches  the  outer  boundary  limit  b*l: 


|(l-exp(-28t) )1/2  «  1 
p 


or  until 


t  "  “28 


Method  (2):  The  method  of  moments  with  j*2. 
Using  j«=2  (a  second  order  approximation) , 


w2(p)  *  *  y 


The  residual  integral  is  then 


/,w2R£dy  * 


/1y{(2y-2y2)b' (t)b(t)  -  2  +  gb2 (t) (l-2y+y2) Jdy  *  0 


The  solution  is 


^WjRgdy  -  |b'(t)b(t)  -  1  ♦  $b2(t)(j|) 


The  differential  equation  is  solved  easily  to  give  the  new 


approximation  for  b(t) 
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CB  *  f(v)  *  |l  “  p 9  n  \/2  I  ^ 

I  [|^Cl-exp{-2Bt))]  1/2  I 

Method  (4) :  The  Galerkin  Method  (p)  *  £  (p)  «  l-2p+  p2 

/‘w.JLdn*/1 U-2p-p2) I (2p-2p2)b(t)b' (t)-2+Bb2 (t) (l-2p+p2) ]dp«0 

•  3  *  • 


Integration  and  solution  of  the  resulting  differential  equation 
yields 


cB  -  f(p) 


1- 


_ x _ 

^|(  1-exp  (-4  8t) 


>] 


2 


(59) 


Method  (5):  Integral  method,  j=l,  but  choose  as  an  approximation 
polynomial  the  function 


cfi  *  f(y)  *  1-sin  ^ 


(60) 


instead  of  equation  (34) . 

Thus 

2 

f’(p)  *  -J  cos  Jp  and  f"(p)  «  j-  sin  Jp 
so  that  since 


(61) 


Rg  -  -pf '  (p)b(t)b*’1(t)-f "  (p)b~2  (t)  +  Sf(p> 


(62) 


we  have 


Rj.— p(-|cosjp)b(t)b“1(t)-(}-sinjp)b"2(t)+B(l-sin|p) 


(58) 


(63) 
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For  comparison  of  the  exact  solution  with  the  approximations,  we 
choose  k  «  100  s T  «  10~2  s,  D  *  10’5  cm2s-1  and  (A#l  *  10~3  moles  l”1. 

We  find  that  if  we  choose,  under  these  conditions,  the  boundary 
X*L  l  3  x  10'3  cm,  that  there  are  no  diffusion  effects  at  the  boundary 
under  these  conditions,  therefore  the  approximate  solutions  are  valid 
using  I*  *  3  *  lo'3  cm.  This  leads  (using  the  relations  in  (19))  to 
the  following  values  of  the  dimensionless  parameters! 


t  -  0.0111 


6  *  90 


(68 
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i'abie  1  lists  the  results  of  each  of  the  approximation  methods  and 
the  exact  solution  at  various  distances  from  the  electrode-solution 
interface.  The  accuracy  of  the  collocation  method  3b  is  impressive  for 
regions  close  to  the  electrode  (or  at  short  times) .  This  method  is 
particularly  interesting  because  of  the  low  order  (j»l)  of  polynomials 
used  for  the  basis  functions.  It  will  be  seen  in  the  next  section  that 
the  use  of  the  more  sophisticated  orthogonal  collocation  method  leads 
to  continuation  of  rapid  convergence  as  the  approximation  order  is 
increased.  The  approximation  is  thus  capable  of  giving  highly  accurate 
results  with  little  computational  effort. 

III.  ORTHOGONAL  COLLOCATION 

A .  Introduction 

The  method  of  orthogonal  collocation  has  as  its  basis  the  existence 

of  what  is  known  as  interpolating  polynomials.  If  we  have  a  real 

function  0(x),  and  we  choose  n  points  x, ,x,- — x_  in  the  interval  [x, ,x  _] 

x  4  n  in 

then  it  may  be  shown  [10]  that  there  exists  some  polynomial  of  degree 
n-1  which  equates  to  0 (x)  at  each  of  the  points  xn.  If  the  number  of 
points  xn  is  increased,  a  new  polynomial  may  thus  -be  found  which  de¬ 
scribes  0(x)  at  these  new  points.  Therefore,  in  the  limit,  0(x)  can 
be  described  at  every  point.  Certain  techniques  may  be  used,  however,  to 
use  a  relatively  small  number  of  xft  (interpolation  points)  and  still 
describe  0 (x)  accurately  at  all  points  in  the  interval.  This  is  most 
readily  accomplished  by  choosing  the  proper  type  of  polynomial.  The 
interpolation  points  are  defined  by  the  polynomial  chosen,  and  are  usually 
the  roots  of  that  polynomial. 

The  purpose  of  orthogonal  collocation  solutions  to  partial  dif¬ 
ferential  equations  in  to  supply  time  dependent  coefficients  to  the 


Table  1.  Comparison  of  results  of  the  several  MWR  for  the  catalytic 

n  *j 

mechanism.  Parameters:  T  *  10~  s,  (A°l  •  10”  M,  k  *  100  s 
D  *  10"5  cm3s”*,  8  *  90. 

I  Concentration,  species  B,  M/10**3 


2/cm 

EXACT 

1 

2 

3a 

3b 

4 

5 

2xl0"6 

.9934 

.9944 

.9954 

.9944 

.9936 

.9930 

.9949 

4 

.9868 

.9889 

.9908 

.9888 

.9872 

.9861 

.9897 

6 

.9802 

.9834 

.9863 

.9832 

.9808 

.9791 

.9846 

8 

.9738 

.9779 

.9817 

.9776 

.9745 

.9722 

.9794 

lxl0~5 

.9673 

.9724 

.9772 

.9721 

.9682 

.9653 

.9743 

2 

.9355 

.9452 

.9546 

.9445 

.9369 

.9313 

.9486 

4 

.8749 

.8920 

.9103 

.8906 

.8758 

.8650 

.8974 

6 

.8177 

.8403 

.8670 

.8383 

.8169 

.8012 

.8464 

8 

.7639 

.7902 

.8247 

.7876 

.7699 

.7398 

.7959 

lxl0"4 

.7133 

.7416 

.7835 

.7385 

.7015 

.6809 

.7459 

2 

.4930 

.5218 

.5934 

.5165 

.4615 

.4230 

.5085 

4 

.2343 

.1977 

.2924 

.1913 

.1287 

.0906 

.2118 

19 


various  powers  of  x  in  the  polynomial,  thereby  providing  0(x)  at  the 
interpolation  points  xn.  In  the  method  described  herein,  the  polynomials 
chosen  are  orthogonal  shifted  Jacobi  polynomials,  and  the  interpolation 
or  collocation  points  are  the  roots  of  those  polynomials.  These 
choices  have  been  shown  [11]  to  be  sufficient  to  describe  solutions  to 
the  types  of  systems  of  equations  that  appear  in  electrochemical 
diffusion-kinetic  problems. 

The  advantages  over  finite  difference  methods  are: 

(1)  Inherent  stability  of  solutions. 

(2)  Greatly  increased  savings  in  computational  effort. 

(3)  Ease  in  changing  programs  from  one  mechanism  to  another. 

A  given  polynomial  of  degree  n-1,  i.e.  Pn_^(x),  may  be  found  if  it 
is  possible  to  solve  a  set  of  simultaneous  equations  in  the  xn  inter¬ 
polation  points,  i.e. 


°n-lxl~1  +  an-2xl~2  +  -  +  °lxl  +  °0  *  ^xi> 

°n-lX2  1  +  an-2x2  2  +  +  alX2  +  “0  “  ® !x2) 


“n-l’C1  +  “n-2xr2  + - +  “lxn  +  «0  *  *(xn> 


If  the  determinate  in  the  x^  of  the  above  system  is  nonzero,  then  a 
unique  solution  exists.  Since  the  xR  are  known,  and  the  orthogonal 
colligation  technique  will  supply  the  o^»  the  solution  of  the  0(xn) 
is  straightforward. 

It  was  mentioned  above  that  as  the  number  of  interpolation  points 

are  increased  in  the  interval,  the  more  accurately  the  polynomial 

interpolates  to  0(x)  at  every  point.  The  exact  error  bound  at  a  point 

jk_,  where  x_yx„  is  given  by  the  relation 
m  m  n  - 
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Error  [0;m]  s  ~  j  Yn(m)  j  2R  (70) 

where  Y  (x)  «  (x-x.)(x-x,) - (x-x  ) ,  and  2.  is  defined  if  we  know  that 

n  id  n  n 

|0(n><x)  |  <  2n  for  all  x  in  lxi*xnK  This  assumes  that  0(x)  has  con¬ 
tinuous  derivatives  0*n*  in  the  interval.  It  is  obvious  that  as  the 
number  n  of  interpolation  points  is  increased,  there  is  a  rapid  decrease 
in  the  error  bounds. 

The  individual  reacting  species'  spectrochemical  absorbances  and 
the  total  charge  consumed  during  the  reaction  may  be  obtained  by  the 
time  integration  of  the  appropriate  concentration  profile  and  the  current, 
respectively.  For  polynomials,  a  class  of  formulas  called  quadrature 
formulas  [12]  exist  such  that  their  solutions  yield  exact  values  of 
the  integral  of  the  polynomial  over  the  designated  interval.  With 
respect  to  a  specified  weight  function  (w(x),  the  formula 

/  n  w(x)0 (x)dx  «  I  Q40(x.  )  (71) 

x1  i«l  1  1 

represents  an  exact  quadrature  formula  for  the  integral  of  0(x)  with 
respect  to  the  weight  function  w(x)  over  the  interval  lxj/xnJ  •  The 
Qj.  are  constants  that  are  determined  once  again  by  methods  dependent 
on  the  type  of  polynomial  that  has  been  used  to  simulate  0(x)  in  the 
interval  lx2»xnl»  Certain  of  these  quadrature  formulas  lead  to  the 
well  known  Newton-Cotes,  trapezoidal,  and  Simpsons'  rule  integration 
formulas.  We  will  be  concerned  here  with  the  highly  accurate  Gauss- 
Jacobi,  Lobatto,  and  Kadau  types  of  quadrature  formulas. 

B.  Discretization  of  the  Differential  Equations 

We  are  interested  in  solving  an  equation  describing  electrochemical 
diffusion  phenomenon  such  as  the  type 
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where 


and 


0. 


(72) 


(73) 


(74) 


He  assume  that  there  is  an  interpolation  polynomial  that  will 
describe  0(x,t)  over  the  interval  (0,LJ.  The  constants  o  in  the  poly¬ 
nomial  described  in  the  previous  section ,  are  now  functions  of  time 
(a(t))  since  we  are  dealing  with  a  function  0  of  x  and  t. 

As  mentioned,  the  orthogonal  polynomial  that  we  shall  use  to  ap¬ 
proximate  0(x,t)  is  the  shifted  Jacobi  polynomial.  These  polynomials, 
Pn(y,6)(x)  are  described  by  the  orthogonality  relation 


/1x6(1-x)yPm(y'6)  (x)  Pm(y,6)(x)  dx  »  0  (75) 

where  x^(l-x)Y  is  the  weight  function;  y  *nd  5  >  **1*  Letting  y  *  6  *  0, 
we  have  polynomials  defined  that  are  suitable  for  the  linear  diffusion 
approximation.  These  are  known  as  the  Legendre  polynomials.  Since 
the  interpolation  polynomial  will  have  order  (n-1)  where  n  is  the  number 
of  interpolation  points  to  be  used  for  calculation  of  the  values  of 
0(x,t),  the  highest  power  of  x  in  the  approximation  will  be  (n-1). 

The  trial  function  for  0  may  thus  be  written: 


n+2  H-l 

0T  -  I  «4(t)x3  1  (76) 

T  j«l  3 

The  n+2  terms  come  from  the  fact  that  x  ■  0  and  x  ■  1  are  also  roots  of 
the  polynomial.  In  terms  of  the  interpolation  (collocation)  points, 
aquation  (76)  becomes : 
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0(x.,t)  «  I  a,(t)xj  ,  i  «  1,2, - n+ 2 

1  j-1  3  1 


This  represents  a  set  of  n+2  simultaneous  equations  in  the  unknown  0 
and  a.  Since  the  Oj  represent  the  elements  of  an  n+2  vector  of  the  time 
dependent  coefficients,  the  0(x^,t)  represent  the  elements  of  a  vector 
in  the  0(t),  and  the  x]|”*  are  the  elements  of  an  (n+2, n+2)  matrix. 
Equation  (77)  may  be  represented  in  matrix/vector  notation  after 
Whiting  and  Carr  [5]  as: 


0(t)  *  a*Q 


at  fixed  x^,  where  the 


We  now  differentiate  (77)  v.'th  respect  to  x  to  get  the  terms  in 
the  original  equation  (76) : 


d  dx^“* 

axW(x't>>  -  ^  •jew  i— 


or  in  general. 


Z-rpriHx,  t)) 
dx'K' 


l  o4  (t) 


a<IOxj-l 


Specifically  then,  at  the  collocation  points, 


a(X>  .  nt2  „ 

SW _  j.i  i 


n+2  .OO-J-l 

I  «,(«  S-rer - 

i-1  3  dx'K' 


,  i  ■  1,2,™ n+2 


In  matrix/vector  notation,  for  k  ■  1  and  k  *  2,  this  is  equivalent  to 
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dg(t) 

dx 


a (t) *C 


(83) 


and 


—  ■■■■  «  a  ( t )  *D 

dx2 


where 


and 


dx?'1 
1 


'ij 


'ij 


a2**'1 

~zr 


(84) 


(85) 


(86) 


The  desired  quantities  are  the  (t) #  and  follow  from  matrix  algebra: 
From  (78)  we  have 


o(t)  -  Q"1  0(t) 
such  that 


dfl(t) 

dx 


d20(t) 

TT 


and 


or 


(87) 

(88) 

(89) 


-  A  0(t)  ,  and 


(90) 


-  B  *(t) 
dx2 

where  A  *  C  and  B  -  DQ~* 

For  a  single  collocation  point#  wa  have 


(91) 
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d0(t) 


n+2 

I  A.  .0  (x .  ,t) 
j*l  x3  3 


(92) 


d20(t) 

dx2 


x. 

l 


n+2 

«  l  B. .0 (x, ,t) 
j=l  13  3 


Equation  (93)  may  be  used  immediately  in  the  diffusion  equation  0 


to  yield 


d0 

dt 


x=x. 


n+2 

I  B. .0 (x . ,t) 
j=l  13  3 


(93) 

*xx 

(94) 


The  diffusion  equation  is  reduced  to  n+2  simultaneous  first  order  dif¬ 
ferential  equations  in  n+2  unknowns,  0(x^,t). 

C.  Integration 


There  are  many  standard  methods  that  may  be  used  for  the  integration 
of  these  simultaneous  first  order  differential  equations.  We  have 
found  113]  that  the  method  of  Calliaud  and  Padmanabhan  [14]  is  probably 
the  fastest  and  most  accurate  for  diffusion/kinetic  equations.  This 
method,  known  as  ISI3,  has  been  shown  [15]  to  be  highly  effective  for 
the  integration  of  stiff  coupled  differential  equations.  The  reader  is 
referred  to  the  original  literature  for  details  of  the  derivation;  the 
pertinent  features  are  given  here. 

In  a  Runge-Kutta  method  of  integration,  many  derivative  approaches 
may  be  taken.  In  the  semi-implicit  ISI3  variation,  the  solution  y  at 
the  n+1  point  to  a  general  differential  equation 


is  given  by 

*n+l  *  *n  +  Vi 


(95) 


(96) 
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Here,  B  is  the  Jacobian  of  f(y)  at  yR,  h  is  the  integration  step  size, 
are  constants,  and  the  I  are  elements  of  the  identity  matrix.  The 
R^  and  c^  are  found  by  comparing  Taylor/power  series  expansions  of 
equation  (95)  with  the  Taylor  series  expansion  of  y  Using  this 

procedure  with  the  constraint  of  making  the  characteristic  root  at  (-») 
equal  to  zero  (a  necessary  condition  for  accurate  integration  of  stiff 
equations),  the  values  in  Table  2  for  the  c^  and  are  obtained. 


Table  2 


c1  0.43586659 

c2  0.75 

e3  -0.27468397 

c4  -0.10562709 

Rx  16/27 

R2  11/27 

*3  1. 

Use  of  these  constants  in  the  above  equations  gives  very  accurate 
integrations  of  differential  equations  possessing  large  differences  in  sub¬ 
sequent  eigenvalues  ("stiff"  equations).  This  is  the  primary  method  used  i i 
this  work,  and  the  Fortran  program  STIFF 3  [15]  is  given  in  the  appendix. 
Included  also  are  programs  for  calculating  the  roots  of  an  orthogonal  poly¬ 
nomial  (collocation  points) ,  the  A  and  B  coefficients  for  the  descretized 
differential  equations,  and  programs  to  calculate  the  quadrature  coef¬ 
ficients  and  perform  polynomial  integration.  The  examples  there  demonstrat 
the  use  of  the  programs. 
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D.  Spline  Collocation 

When  dealing  with  equations  containing  large  homogeneous  rate  con¬ 
stants  and/or  trying  to  solve  equations  at  very  short  times,  it  is  obvious 
that  the  use  of  a  low  collocation  approximation  order  will  not  be  sufficient 
to  simulate  the  response  accurately.  The  rapidly  changing  response  simply 
might  take  place  before  the  first  interior  interpolation  point.  This 
is  the  same  problem  that  presents  itself  in  finite  difference  schemes 
and  to  which  Joslin  and  Pletcher  [2]  have  treated  for  that  technique. 

In  collocation  techniques,  there  are  several  good  ways  to  eliminate  this 
problem  while  maintaining  the  use  of  only  a  few  collocation  points  for 
the  simulation.  One  such  method  is  the  automatic  choice  of  the  pre-sum¬ 
mation  factor  as  discus  ->d  in  section  (F.b.).  Here,  another  simple  option 
is  discussed:  the  method  of  global  splines. 
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There  are  techniques  for  dividing  the  interval  of  interest  into 
many  subintervals,  and  then  performing  a  standard  collocation  technique 
on  each  subinterval  [15] .  The  solutions  at  the  outer  boundary  of  each 
subinterval  become  the  internal  boundary  conditions  for  the  next  sub¬ 
interval.  The  procedure  continues  until  the  original  outer  boundary 
(or  last  subinterval)  has  been  treated.  The  number  of  equations  is 
raised  by  a  factor  equal  to  the  number  of  subintervals.  This  technique 
then  becomes  similar  to  simply  raising  the  number  of  mesh  points  in  a 
finite  difference  scheme,  although  one  finds  that  for  similar  accuracy 
between  the  two  methods,  this  spline  technique  is  still  significantly 
faster. 

The  diffusion  boundary  concept  provides  a  very  simple  but  highly 
efficient  spline  method.  Initially,  the  original  spacing  between  the 
collocation  points  is  compressed  so  that  their  total  spanned  interval  is 
just  slightly  larger  than  the  region  where  the  profile  is  changing  most 
rapidly.  As  the  width  of  the  change  increases,  the  concentration  value 
is  tested  at  the  last  internal  collocation  points  (N  or  N+l) .  As  the 
concentration  there  begins  to  become  larger  (or  smaller  if  appropriate) 
than  the  boundary  condition  for  that  species,  then  the  distance  between 
the  collocation  points  is  expanded,  and  the  procedure  repeated.  In  this 
new  repetition,  solutions  at  all  times  less  than  the  time  at  which  the 
expansion  was  made  are  discarded  as  inaccurate.  The  procedure  is  con¬ 
tinued  until  the  final  desired  time  has  been  reached.  Solutions  of 
exceptional  accuracy  are  obtained.  For  instance,  the  maximum,  relative 
error  in  a  catalytic  mechanism  simulation  under  a  chronoamperometric 
experiment  where  the  rate  constant  was  100  s“*  and  the  integration  was 

.  .n 

from  10  s  to  20  s,  the  maximum  relative  error  was  0.07  percent. 
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The  procedure  for  implementing  this  "spline  point”  method  is  quite 
simple  to  execute.  The  subinterval  is  normalized  to  a  (0,1)  interval 
by  a  variable  transformation,  end  the  normal  integration  procedure  is 
followed.  However,  the  value  of  the  concentration  at  the  outer  points 
is  tested  continuously  as  mentioned.  The  test  value  is  dependent  on 
several  factors,  including  the  rate  constant,  time  of  integration, 
number  of  species,  etc.  It  was  found  empirically  that  in  using  6  internal 
collocation  points  (6th  order  Legendre  polynomial)  for  the  ECE/DISP 
mechanism  under  a  chronoamperometric  experiment,  that  a  branch  test  at 
the  N+l  point  of  0.001  concentration  units  was  sufficient  to  maintain 
high  accuracy  even  at  very  high  rate  constants  (up  to  106  first  order 
and  10*®  second  order)  when  compared  to  finite  difference  schemes  taking 

4 

computing  times  up  to  10  times  as  long  and  extremely  fine  mesh  sizes 
(high  memory  usage) . 

Consider  the  general  partial  differential  equation 


(100) 


in  the  interval  (0,1).  If  we  desire  to  insert  a  spline  point  at  some  value 
x,  such  that  0  <  x#  <1,  then  we  would  simply  make  the  variable  trans¬ 
formation 


(101) 


in  our  equation  to  renormalize  the  boundaries.  The  equation  then  reads 


«£.  *  O 
ft  x2  fij2 


(102) 
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Discretisation*  then*  leads  to 


60 

6t 


s. 

1 


— y  I  B-.|F(x,  *t) 

V  3'1  3 


(103) 


Integration  of  the  set  of  equations  (103)  is  performed  as  usual.  When 
the  branch  is  made  to  a  wider  interval,  xg  is  changed  and  the  process 
repeated. 

Spline  techniques  are  given  in  example  3  in  section  C. 


E.  Discretization  Examples 

Example  1:  The  reversible  charge  transfer  mechanism  under  a  potential 

step  to  a  region  where  the  kinetics  are  diffusion  con¬ 
trolled. 

ne~ 

A  B  (104 

The  descriptive  equations  and  boundary  conditions  are: 


«T 


Mil  .  B  «!l|l 

4T  B  {X2 


<a1o.t  •  -  IB1x,o  ■  0 


lA1x,o  ■  <a>.,t  -  <A*J 


Inserting  the  dimensionless  parameters 

L5  E  ,  and  cA  -  » 


(105 


(10€ 

(107 

I 


!10' 


we  have 
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After  partial  expansion,  we  obtain 


dcB  N+l 

Ht~  *  Bi,lcB<0,t*  +  Bi,N+2CB(1,t)  +  Bi,jCB(xj#t) 


(110) 


(111) 


cA(0,t)  *  cB(l,t)  •  cfi(x,C)  *  0 


(112) 


cA(x,0)  «  cA(l,t)  *  1 


,CA 

X1 


xi*o 


Prom  (93) ,  we  then  have 


(113) 


(114) 


N+2 

£  B.  ,c  (x,.t> 

j«l  A  3 


(115) 


N+2 

I  B.  .c  (x.,t) 
+*i  +»J  »  3 


(116) 


Bi,lCJi(0,t)  4  Bi,N+2CA(1,t>  +  ^l2  Bi,jCA<Xj't) 


(in: 


The  known  boundary  conditions  (112)  and  (113)  for  cA(0,t),  cA(l,t)  and 
cB(l,t)  nay  be  entered  at  this  point  to  give 

dcA|  N+l 

an,  *  ®i,N+2  ♦  *  Bi,jcA<*j’t)  \ 


ar  _  *  *  *  Bi,jcB<xj't) 
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Mote  that  we  must  use  some  other  means  to  eliminate  cfi(0,t).  This  may 
be  done  by  discretizing  the  boundary  condition  (114)  with  (92). 


xr° 


Xj*0 


(121) 


i  •  «»  •*  •  •• 

I  A.  *  *  I  A.  •<Cp(Xj  (t) 

A  3  j*l  A»3  B  3 


(122) 


*lAl#lCA(0't)‘i'Al,N+2CA<1,t)+^2Ai,jCA<Xjrt)1  *  Al,lCB(0't) 


+Al,N+2CB(1#t)',‘^2Ai,jCB(xj#t) 


(123) 


Inserting  the  known  boundary  conditions  (112)  and  (113)  again  for  cA(0,t), 


cA(l,t)  and  cfi(l,t)  we  have 


•[A.  M.0+  I  A.  .c_(x.,t)1  «  A.  .cft(0,t)+  I  A.  ,c_(x.,t) 
l,N+2  A»3  A  3  1,1  B  *'3  B  3 


(124) 


Solving  for  cfi(0,t). 


:B(0,t)  -  -  lAi#N+2+^2Ai,j(cA(xj,t)+cB(xj,tn 


(125) 


This  explicit  value  may  now  be  substituted  in  (120)  to  yield 
dcR  B.  ,  N+l 

3t"  *  "  A^“[  lAl»N+2+j^2Ai#j(CA(Xj't>+CB(Xj,t)1 


+  I  B,  .c_(x.,t) 
j«2  1,3  B  3 


(126) 


Mow  the  problem  has  been  reduced  to  the  simultaneous  solution  of 
2N  equations,  i.e.  (119)  and  (126)  in  order  to  obtain  c^(t)  and  cp(t ) 
at  each  of  the  collocation  points  x^. 
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We  will  depart  from  the  normal  format  at  this  point  to  completely 
develop  the  program  necessary  to  solve  this  particular  basic  example. 

In  the  later  examples,  only  the  equations  necessary  to  be  programmed 
will  be  given. 

In  this  example,  we  use  a  Legendre  polynomial  of  order  six  (six 
internal  roots  or  collocation  points) .  We  know  the  concentration  at 
both  boundaries  for  the  A  species  at  all  times,  and  we  know  the  con** 
centration  of  B  at  the  outer  boundary  at  all  times.  The  concentration 
of  B  at  the  electrode  surface  may  be  found  easily  if  we  know  all  of 
the  other  A  and B  concentrations  at  each  internal  point  as  one  can  see 
from  (125) .  For  a  large  savings  in  computational  time  and  programming 
ease,  it  will  be  simpler  to  simply  solve  for  the  six  internal  con¬ 
centrations  of  A  and B  (12  equations)  instead  of  at  all  8  points  (16 
equations).  The  procedure  calls  for  the  writing  of  3  short  subroutines, 
FUN,  DFUN,  and  OUT.  In  this  example,  FUN  consists  of  defining  the 
differential  equations  and  returning  their  values  as  F(I),  I  •  1  to  6 
for  the  A  species  (119),  and  F(I+6),  I  ■  1  to  6  for  the  B  species.  We 
will  use  Y  (I) ,  I  «  1  to  6  for  the  six  concentrations  of  species  A  at 
the  internal  points,  and  Y  (1+6) ,  I  *  1  to  6  as  the  six  internal  con¬ 
centrations  of  species  B.  From  (119),  therefore,  we  have 
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the  i  and  j  indices  on  the  right  hand  aide  of  the  equations  on  the  A 
and  B  coefficients  are  simply  changed  to  i  +  1  and  j  +  1.  We  are 
simply  xefering  to 

dcA 

F(I)  *  i  »  2,7  E4 

xi 


as 


F (I)  «  ~  i  *  1,6  E5 

xi 

The  j  index  on  the  summation  thus  is  changed  also  from  j  *  2  to  7  to 
j  *  1  to  6  to  maintain  consistency. 

Partially  transforming  to  FORTRAN  language,  we  have 

6 

F(I)  »  B  (1+1,  8)  +  E  B(I+1,J+1)*Y(J)  £6 

j-1 

To  get  rid  of  the  summation  symbol,  we  use  the  following  loop 

DO  1,  I  «  1,6 
TDDER1  *  0. 

DO  2  J  «  1,6 

2  TDDER1  «  TDDER1  +  B (1+1 , J+l) *Y (J) 

1  F(I)  ■  B (1+1,8)  +  TDDER1 

Thus  the  6  differential  equations  are  defined  in  subroutine  FUN  and 
will  be  returned  to  the  main  program  for  simultaneous  integration. 

Programming  for  the  six  differential  equations  F(I+6),  I  *  1  to  6 
for  the  six  internal  concentrations  Y(X+6)  for  the  B  species  follows 
a  similar  line  of  reasoning.  We  also  "re-index"  for  programming 
simplicity  in  the  same  manner.  (Replace  i  and  j  on  right  hand  side  of 
(126)  on  the  collocation  coefficients  A  and  B  by  i  +  1  and  j+l  and 
change  the  summation  index  from  j  ■  2  to  7  to  j  ■  1  to  6.) 
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F  (7) 


~^IA.  I  A,  ...  (c.  (x.,t)+cR(x,,t)) 
Al,l  1,8  j*l  2»J+1  A  3  B  3 


+AB2,  j+^B^j'1*1 


FU2) 


q^7lAl  ,  8+^1A7 ,  j+l  {cA (x j ' +CB (X j #  t} } 


+  I  B.  ...cR(x.,t)] 
jel  3  0+1  B  3 


Partial  transformation  to  FORTRAN  gives 

6 

F { 1+6 )  *  -B (1+1 ,1) /A (1,1) *  (A  (1, 8)  +  I  A  (1+1 ,  J+l) * (Y (J) +Y (J+6) ) 

J«1 

6 

+  I  B(I+l,J+l)*Y(J+6) 

J=1 


Again  the  summation  signs  may  be  eliminated  by  the  following  loops 

DO  1  I  »  1,6 
TDDER2  *  0. 

TDDER3  «  0. 

DO  2  J  «  1/6 

TDDER2  »  TDDER2+A(I+l,J+l)*(Y(J)+Y(J+6)) 

2  TDDER3  *  TDDER3+B (1+1, J+l) *Y (J+6) 

1  F (1+6)  «  -B (1+1) /A (1/1) * (A (1, 8 ) +TDDER2 ) +TDDER3 


Thus,  the  entire  subroutine  defining  both  sets  of  differential  equations 
for  the  A  and  B  species  would  appear  as  follows: 

(  1)  SUBROUTINE  FUN  (Y,F) 

(  2)  IMPLICIT  REAL* 8 (A-H , 0-F ) 

(Implicitly  invoking  double  precision  mathematics  for  all  real 
variables) 

(  3)  DIMENSION  Y(16),  F(16) 

(  4)  COMMON  A(30/30),  B(30,30),  ROOT(30) 

(The  collocation  coefficients  A  and  B  and  the  polynomial  roots 
ROOT  have  been  defined  by  subroutine  JCOBI  and  DFOPR.  They  are 
passed  to  FUN/  when  needed/  by  the  COMMON  block  statement.) 

(5)  DO  1  I  -  1,6 
(  6)  TDDER1  ■  0. 

(  7)  TDDER2  -  0. 

(  8)  TDDER3  ■  0. 

(9)  DO  2  J  -  1,6 
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(10)  TDDERl  «=  TDDERl+B (1+1 , J+l)  *Y  (J) 

(11)  TDDER2  «  TDDER2+A(I+l,J+l)MY(J)+Y(J+6)) 

(12)  2TDDER3  «  TDDER3+B (1+1 , J+l) *Y (J+6) 

(13)  F (I)  -  B ( 1+1 , 8 ) +TDDER1 

(14)  1  F (1+6)  «  -B (1+1 ,1)/A (1/ 1) * (A (1, 8) +TDDER2) +TDDER3 

(15)  RETURN 

(16)  END 


The  next  subroutine  that  must  be  defined  is  DFUN.  DFUN  is  used 
by  the  integration  subroutine.  DFUN  must  return  the  values  of  the 
derivatives  of  F(I)  and  F(I+6)  at  each  one  of  the  collocation  points 
xi  as  DF (I)  and  DF(I+6),  I  «  1  to  6. 

The  differentiation  is  seen  to  be  straightforward  upon  expansion 
of  one  of  the  F(I)  terms,  say  F(l) 


DF (1,1) 


Ell 


DF (1,2) 


d (F (1) )  .  k 

dcA^2't)  *'3 


E12 
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These  are  programmed  into  DFUN  as  follows 

DO  1  I  «  1,6 
DO  1  J  »  1,6 
1  DF (I , J)  *  B (1+1, J+l) 

The  F(I«6)  are  differentiated  similarly.  For  instance. 


p(7)  «  -  -  ' .  ^Al,8+A2,2  ^A^l'^^B^l'^  *+A2,3^CA^X2't^+CB*X2't*  * 

1  f  1 

+  A2,4(cA(x3,t)tcB(x3#t))  +  A2#5(cA(x4,t)+cB(x4,t))  + 

A2,6(CA<X5't>+CB(X5't)>  +  A2,?(CA(x6,t)+CB(x6't)  +  B2,2CB(Xl't) 
+  ®2,3CB*X2't’  +  B2,4Cb'X3'4'  +  B2,5CB*X4't*  +  B2»6CB*XS'4* 

+  B2,7°B(Vt)]  E1 

"e  desire  ■  i-1-6  E1 

*  3 

Note  that  the  cfi(Xj,t)  are  the  Y(I+6),  so  that  we  use  the  nomenclature 
DF (1+6 ,  J+6)  for  the  B  species. 

DFI7‘7>  -  dCgix^.t)  *  -  £^7*2,2  4  ®2,2  El1 


DF(7'8)  *  icglxJTtT  *  •  A^7  A2,3  +  B2,3 


OPf7  nr- .  jumu .  A  +  B 

OF{7'12/  *  dcB(Xg,ti  Kltl  A2,7  +  B2,7 


or  in  general, 

DF (1+6, J+6)  -  dc~  IxT^y  "  "  "A^TTty^  *A(I+l,J+l)+B(X+l.tJ+l) 

®  j 

This  is  programmed  exactly  the  same  as  the  DF(I,J)  above.  The  entire 
subroutine  then  appears  as  follows 
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(  1)  SUBROUTINE  DFUN  (Y,DF) 

(  2)  IMPLICIT  REAL  *  8  (A-H,0-Z) 

{  3)  DIMENSION  DF(30,30) 

(  4)  COMMON  A (30, 30),  B(30,30),  ROOT(30) 

(5)  DO  1  I  «  1,6 

(6)  DO  1  3  *  1,6 

(  7)  DF (I ,  J)  *  B(I+1,J+1) 

(  8)  1  DF (1+6, J+6)  «  -B (1+1 , 1) /A (1, 1) *A (1+1, J+l) +B (1+1 , J+l) 

(  9)  RETURN 

(10)  END 


The  OUT  subroutine  is  called  at  predetermined  time  increments  from 
the  integration  subroutine  STIFF3.  When  OUT  is  called,  the  0  values 
for  the  concentration  profile  (Y(I)  and  Y(I+6),  I  *  1  to  6)  are  available 
for  output.  Four  parameters  are  transferred: 

X  ■  the  current  value  of  the  time 

Y  *  the  current  value  of  the  concentration  at  each  collocation  point 
IH  *  number  of  bisections  that  occurred  before  successful  integration 
Q  «  stepwise  acceleration  integration  factor. 

In  general,  IH  and  Q  are  of  little  interest  for  our  purposes. 

One  thing  that  may  be  done  first  in  OUT  is  the  calculation  of  the 
unknown  boundary  condition— i.e.  the  value  of  the  B  species  at  the 
electrode  surface  cfi(0,t).  For  convenience  we  call  this  value 
Y(13).  From  (125),  after  re-indexing, 

•1 


*(13)  .  cB(0,t).-~  t^,8+J  *L+14+ltVYt>+CB(Vt)) 


E19 


We  perform  this  calculation  in  OUT  since  all  of  the  cA  and  cfi  (Y(I)  and 

Y(I+6))  are  available.  Partial  conversion  to  FORTRAN  gives 

DO  1  I  •  1,6 
TDDER4  ■  0. 

DO  2  J  -  1,6 

2  TDDER4  -  TDDER4+A(I+1, J+l) * (Y (J)+Y (J+6) 

1  Y (13)  ■  -  i/A (1,1)* (A (1,8) _TDDER4 ) 

The  rest  of  the  OUT  subroutine  may  be  used  for  calculating  the  current 

or  spectrochemical  absorbance,  etc.,  as  well  as  defining  the  output 

format. 
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For  instance,  the  flux  at  the  electrode  surface  is  given  by 


Partial  expansion  gives 


dx  0  *  Al,lCA(0't)+Al,8CA(1't)+j*2Al,jCA(Xj't> 


A.  1  A.  jC, (x. ,t) 

*,*>  A  J 


after  insertion  of  cA(0,t)  *  0  and  cA(l,t)  *  1. 

Switching  to  the  new  index  for  programming ,  we  have 


FLUX  «  A(l, 8)  4  I  A(l,J4l)*Y(J) 

J«1 

Eliminating  the  summation  sign,  we  have 

TDX  ■  0. 

DO  1  J  *  1,6 

1  TDX  «  TDX4A(1,J41)*Y (J) 

FLUX  «  A ( 1 , 8 ) 4TDX 

So  to  calculate  the  surface  concentration  of  the  B  species  at  all  times, 

and  to  calculate  and  print  the  flux  as  well  as  the  profiles  for  each 

species,  the  subroutine  would  appear  as  follows: 

(  1)  SUBROUTINE  OUT (X,Y,IH,Q) 

(  2)  IMPLICIT  REAL* 8 (A-H,0-Z) 

(  3)  DIMENSION  Y (30) 

(  4)  COMMON  A(30,30) ,B(30,30) ,ROOT(30) 

(5)  DO  1  I  «  1,6 

(  6)  TDX-0. 

(  7)  TDDER4-0. 

(8)  DO  2  J  -1,6 

(  9)  TDDER4«TDDER44A(I41,J41)*(Y(J)4Y(J46)) 

(10)  2  TDX«TDX4A(l,J4l)*Y(J) 

(11)  Y(13)  —  l./A(l,l)*(A(l,8)4TDDER4) 

(12)  1  FLUX-A (1,8) 4TDX 

(13)  WRITE (6, 100 )X, FLUX, (Y(I) ,1-1,6) ,Y(13) , (Y(I46) ,1-1,6) 

(14)  100  FORMAT ( '  TIME-JE20.15,/, *  FLUX- ' ,E20 .15,/, 4X, 

1  * 0.000^0000000 J 3 (F16.ll, IX) , /IX, 3 (F16.ll, IX) , 3X, 

1  *  1. 00000000000 J//,4X, F16.ll, IX, 3 (F16.ll, IX) ,/, 

1  IX, 3(F16. 11, IX) ,3x,*0. 00000000000’) 

(15)  RETURN 

(16)  END 


The  main  program  for  the  solution  of  these  equations  would  appear,  for 
example,  as  follows: 


IMPLICIT  REAL*8 (A-H  O—S) 

DIMENSION  DIF1 (30) ,DIF2 (30) ,DIF3 (30) ,VECT{30) ,F(30) ,FOLD(30) 

DIMENSION  Y0LD1 (30) ,YA (30) ,YK1(30) ,YK2(30) ,YK3(30) 

DIMENSION  DF (30,30) ,DFOLD(30,30) ,W(30) ,Y(30) ,YOLD(30) 

(The  above  defines  the  subroutine  variable  arrays.) 

COMMON  A(30,30) ,B(3Q,30) ,ROOT(30) 

EXTERNAL  FUN,DFUN,OUT 
N«6 

(Defines  order  of  polynomial  used  and  hence  the  number  of 
internal  collocation  roots.) 

N0=1 

Nl**l 

(Specifies  that  the  X-0  and  X-l  boundaries  will  also  be  roots, 
(see  JCOB1 ) ) 

AL*0 

BE=0 

(Specifies  that  a*B=0  and  hence  an  orthogonal  Legendre  poly- 
nomial  will  be  used  (see  JCOB1).) 

ND=30 

(Indicates  maximum  vector  dimension  size  for  this  program  -  note 
this  is  the  array  sizes  in  the  DIMENSION  statements.) 

CALL  JCOBI (ND , N , N0 , NL , AL , BE , DIFl , DIF2 , DIF3 , ROOT ) 

(Calculates  the  roots  of  the  Legendre  polynomial,  i.e.  the  X..) 
DO  4  I  «  1,8  * 

4  A(I,J)«VECT(J) 

(Calculates  the  A  matrix  collocation  coefficients  and  places 
them  in  the  A(I,J)  matrix.) 

DO  5  I  «  1,8 

CALL  DFOPR (ND ,N ,N0 ,N1 , I , 1 ,DIF1 ,DIF2 ,DIF3 , ROOT , VECT) 

5  DO  5  J  «  1,8 
B(I,J)«*BECT(J) 

(Calculates  the  B  matrix  collocation  coefficients  and  places 
them  in  the  B(I,J)  matrix.) 

X0*O.D0 

Xl*2O.D0 

(Sets  the  starting  and  ending  times  for  integration.) 

EPS-l.D-06  , 

(Sets  integration  error  limit  at  10  .) 

DO  6  I  -  1,6 

Y(I)«1.D0 

Y(I+6)«0.D0 

(Sets  initial  conditions  for  A  and  B  species.) 

W(I)«1.D0 

6  W(I+6)  *  1.D0 

(Sets  weight  factors  for  results  at  1.) 

H0-1.D-O6  _g 

(Sets  initial  value  for  integration  interval  time  step  at  10 
units.) 

CALLSTIFF 3 ( 12 , 30 , 10 , FUN , DFUN , OUT , X0 , XI , H0 , EPS , 
1W,Y,Y0LD,Y0LD1,IP,YA,YK1,YK2,YK3,DF,DF0LD,F,F0LD) 

(Begins  integration  of  the  equations  contained  in  FUN.  There 
are  12  equations  to  be  integrated  and  OUT  will  be  called  every 
10  integration  steps.) 

STOP 

END 
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These  four  user  defined  programs  (MAIN,  FUN,  DFUN,  OUT)  are  then 
compiled  with  the  six  library  subroutines  STIFF3,  JCOB1,  DFOPR,  SIRK, 
LU,  and  BACK  (which  perform  the  integrations). 
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Example  2:  The  ECE/DISP1  mechanism  under  a  potential  step  to  a  region 
where  the  kinetics  are  diffusion  controlled. 


A - *  B 


B  C  +  H 


n2e” 


EJ  «  EJ 


(127) 


B+C  — k  A+D  Homogeneous  e  transfer,  equilibrium  constant  K. 


3 

D  ^  E 
k. 


The  above  mechanism  has  been  investigated  by  Saveant  et  al.  [16]  and 
Bewick  et  al.  [17]. 

The  equations  to  be  solved  are: 


W  *  da  ^ +  VB>  'c> 


-KjWICl 


^  *  °c  "I1  +  kl,B1  *  VB1  IC) 

ox 

^  ■  DD  ^  +  k2*B)  [C1  -  k3IB)  +  k4IEJ 
OX 

^  *  °E  ^  +  kj'D>  -k4'E> 


The  boundary  conditions  are: 


(128) 


33 


IAJ0  T  «  IB]^  T  -  IC).#T  «  c  0 


(133) 


<B1X,0  "  IC1X,0  *  ^X,0  '  lElX,0  *  0 


(134) 


6A 

**  0/T 


(135) 


Using  the  dimensionless  parameters  t  *  — =■, 


X  r 

v  4 


and 


(136) 


k2 

~  IA°) 
K1 


(137) 


(138) 


(139) 


the  equations  (128-132)  and  boundary  conditions  (133-135)  are  converted 
to  the  following  dimensionless  forms: 


T  +  6cbcc 


(140) 


BcB  *  CB 

— £  «  — ®  -  c  -  Be  c 

6 1  6x2  B  B  c 


(141) 


+  CB  *  BcBCC 


(142) 


+  ecnc«  -  YCn  ♦  6c. 


(143: 


34 


6c, 


-c  s2«V 

-  *  °  - 5  ♦  *CD  -  «c 

fit  fix^  D  £ 


(144) 


cA(0,t)  «  cB(l#t)  *  Cc(l,t)  «  cD(l,t)  «  Cc(l,t)  -  0 
Cg(x,0)  *  cc(x,0)  »  cD(x,0)  •  c£ (x,0)  *  C 


(145) 

(146) 


dCA 

.  Sce 

dcC 

.  dcD 

nr 

V  *■ 

x«0  '  ®" 

1 

Xl-o  ’  ^ 

Xl“° 


(147) 


llote  that  the  normal  assumption  of  setting  all  of  the  diffusion  co¬ 
efficients  equal  to  each  other  is  not  necessary  in  orthogonal 
collocation  simulations,  and  is  simply  done  here  for  brevity. 

Precisely  the  same  technique  is  used  as  in  the  first  example  above 
to  discretize  the  equations.  The  results  are: 


dc. 


N+l 


,  -  otni.H+J+^,Bi3CA(Xj't>1  +  ecB(*i.t)cc(xi,t) 

^  J  * 


(148) 


dCB 

nr 


B 


it! 


N+l 


a  I-  IA,  I  A,  .(c.  (x.,t)+c.,(x.,t))] 


X. 

1 


l"l,N+2  Tj*2"l,j'wA'*j,WTWB'‘kj 


N+l 


dcc 

W 


nr 


+  ^z  Bj, jCfi (Xj ,t) ]  -  cB(xi,t)-8cB(xi,t)cc(xi,t) 


11+1 

*  a  I  Bi^cc(x^,t)+cB(x;j,t)-BcB(xi,t)cc(xi,t) 

xA  j“2 


a*?) 


(150) 


Bia  N+1 


N+l 


*  IT2""  l  1  Ai  ^(cr(x,,t)+cn(x.,t))l  ♦  IB,  .cn(x.,t)l 

Al#l  j«2  c  3  °  3  j«2  4»3  u  J 

+  BcB(xi,t)cc(xi,t)  -  ycD(xi,t) +ficc(:c^»t) 


(151) 


$ 

f'r 

A 

fjj 

.  J  V* 

>V,J 

VI 

* 

**jj 

';P 

It 

1 

9  / 


4 


0 


I 
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dCE 

1ST 


Xi 


B,  .  N+l  N+l 

ol-  1tlL  l  I  A  ,C_(x.,t))  ♦  IB.  *c_  (x  •  .)) 
*1,1  j«2  1,3  E  3  j«2  *'3  E  3,t 


♦  YCDt*xi,t)-5cE(xi,t) 


(152) 


This  last  equation  makes  use  of  the  fact  that  the  flux  of  the 
species  £  at  the  electrode  surface  is  zero,  i.e. 


dcE 

sr 


N+2 

l  A,  .C_  (x .  ,t)  *  0 


x.«0  j=l 


1# j  *  j 


N+l 


*l,lCE(0#t)  +  Al,N+2CE(1,t)  +  :jJ2Al,jCE<Xj#t)  *  0 


(153) 


US') 


Since  c£(l,t)  *  0,  we  have 


c£(0,t) 


ll.l 


I  AlfjCE(x^,t) 


(155) 


which  is  substituted  into  equation  (144)  after  initial  discretization. 
Spline  Modified  Collocation  Examples 
Example  3:  Chronoamperometric  single  pulse 

(a)  Reversible  charge  transfer  under  diffusion  controlled  kinetics. 


-ne- 

A 

B 

(156) 

+ne“ 

*  .  TD 

t  *  t 
ir 

j,  ,  X  2»X_  c  m  f  A| 

1*  XB  A  [A  J 

CB  ■  ■fl'T  ‘’■DA  +D8 

(157) 

cA(0,t) 

*  cB(x,0)  «  cD(l,t)  »  0 

(158) 

eAU,t) 

■  cA(x,0)  ■  1 

(159) 

£ 
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xl*0  /xi*° 


7“2  lDi,N+2+.^®i,jCA(zj#t)1 
xs 


dcR  ,  H+l 

ar  *  ^2  IBi.i4J*2Bi.ieB<vt)1 


(b)  Reversible  Charge  Transfer-Reversible  Chemical  Reaction  (EC*  ) 


-ne- 


*-?  -i  * 


2-  K  -  ^  k-V^  B'(i?k)(^ 


»-■{&■  cb"H^T 


D  ‘  da  '  db  -  Dc 


cA(0,t)  «  cfi(l,t)  «  cc(l,t)  «  cB(X,0)  *  cc(X,0)  -  0 


c.  (lft)  *  c.(XfO)  *  c. 


Xj*0 


,5X  h-° 


.  N+l 

7?  lBifK+2+j^BifjCA(Zj't)1 


dcB  1  (  1  N+1 

ar  ■  7?  j“  j^“IAi#iH2+ji2Ai,jlcA<8jft,+cB(,j't,n 

N+l  ) 

+  ^2Bi,jcB<zj#t)  I  "  6IcB(2i,t)-Kcc(2i#t)] 


(160) 


(161) 


(162) 


(163) 


(164) 


(165) 


(166)  f 


(167) 


(168) 


(169) 


*- 
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dc.  ,  (  B.  .  N+l  W+l  ) 

IT  *  -*5  r  ir^  *  A1  (2.rt)  +  I  B,  (2.,t)> 
dt  2i  x  2  (  *1.1  j«2  C  3  j-2  1'3  C  3  J 


+  0 ICgtz^ft)  -  Kcc(2i#t>l 


(170) 


(c)  Reversible  Electron  Transfer  -  Dimerization 


D  2B  ^  C 


(171) 


*-! 


v  9  ■£,  v 

h  kf[A°)  * 


(A*)kf  +  kfc 


C.  *  j^lr  c„  *  tv-1'.  c_  *  •  '9 1-.-  D  *  D.  *  D„  *  D„ 

A  [A  j  B  (A°J  C  (A®T  A  B  C 


0^(0, t)  «  cB(l#t)  *  Cc(l,t)  *  cfi(x,0)  ■  Cc(xr0)  *  0 


c« (1 » t)  *  c.(x,0)  «  C- o 


(172) 


(173) 


(174) 


(175) 


Xj*0 


Xl*° 


(176) 


dcB  i  (  ®i  i  r  H+1  1 

ar  ^  {-  s^7  ,cA<*j't)+eB(,j*t»J 

•HI  j  , 

+  jcB(Zjft)j  -  0IcB  (zi,t)-Kcc(si,t)J 


B,  ,  11+ 1 


uwp  %  vj  «  «Ti  nTi, 

jrr-  »  — y  1“  g  Y  E  A,  jCr(zH,t)  +  I  B,  Hcc(z.,t)) 

*i  V  *1'1  3*2  1,3  C  3  j-2  1,3  C  3 


+  Pt5CB2(zift)-Kcc(zA,t)J 


(177) 


(178) 


(179! 
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Reversible  Charge  Transfer  -  Disproportionation 


+ne  Kf 

- — *  b  a  +  b  ;p  C 


TD  X  .  X  y  ~  b 

t  *  x  “  l  2  xgK  cA#kf 


k«tA*]kf+kb  B 


(i«)(d  ) 


5a  “  TIFT  cb  *  Cc  *  *^T 


c.  (0, t)  -  c  (l,t)  *  cc(l,t)  *  CB(1,0)  *  CC(1'0)  *  0 


C.(lrt)  -  CA(Z,0)  -  IA#1 


L«0  "  Vs  /  V°  VX  /Xl“° 


sr|  *  "75  1Bi.lH2+I'ijBi.jCA(Sj’tn  ‘  6leAISi,t,CB<Si't> 


zi  *s 


*ScC*2i't)5 


‘i  *‘s 

11+1 
+  ID 
j*2 


A  {-  ^7  [Al.N+2+”i12Al,J1CA,2i't)+CB(23'tn] 
-*s 

;1®i#  jCB(Zj#t)  |  “ 


B.  ,  11+1 


3T  7  -  7?  l“  A^I  ^2Ai,icc(Vt)+  jiaBi.JCc<,i,t>1 


+  8lcA(si,t)cB(si,t)-Kcc(zi#t)] 


(ISC) 


(181) 


(102) 


(183) 


(104) 


(105) 


(106) 


(107) 


(183) 


(e)  Reversible  Charge  Transfer  -  Catalytic  Reaction 

(The  FORTPAll  listings  for  this  exanple  are  given  in  the  Appendix.) 
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-ne~  k 

" - -  B  A  +  C 

(189) 

+ne”  IK) 

-r  x * e 

kL2 

b  n  ^  T\  y  H  y  TJ 

D  A  B  C 

(190) 

\  *  {a^r  CB  ■  TaT 

cc'M 

(191) 

cA(0,t)  *  cB(l,t)  *  cfi(z,0)  «  0  (192) 

cA(l,t)  «  cA(z,0)  *  cfl(0,t)  *  cA(zi,t)  +  Cg (z^  *t)  «  1  (193) 


(194) 


(195) 


(196) 


(197) 


(193) 


where  ie  the  Kronecker  delta  (■!  ifi«j;«Oifi^j) 
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Example  4:  Spherical  Diffusion  -  Chronoamperometric  Single  Pulse 

(a)  Reversible  Charge  Transfer 
-ne~ 

A  B  (199) 

4ne 


The  equations  for  the  system  are: 


(200) 


tA] 


VT 


(201) 


W-.T  *  lAWn,0  “  ^ 


(202) 


6  [A] 
6R 


R«R, 


'  6[B) 
6.R 


R»R, 


(203) 


where  R  is  the  radial  distance  from  the  center  of  the  spherical  elec¬ 
trode,  and  Rq  is  the  radius  of  the  electrode.  We  choose  the  following 
convenient  variable  transformations: 


t 


12 _  r  -  ^  D 

<L-R0)2  l'R0 


Effecting  these  transformations  on  equation  (200)  we  obtain 

ffA  .  +  2(L-R0)2  fA 

fit  fir2  r(L-R0)2+R0(L-R0)  fir 


(204) 


(205) 


or,  after  simplifying, 


-  + 

2 

4ca 

6t 

17 

r+8* 

6r 

(206) 


where 
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e* 


T 


Also,  we  have 


cA(0,t)  *  0 
cA(l#t)  -1 


(207) 


(208) 

(209) 


6cA 

Tr~ 


rl 


6cB 

"  IT 


(210) 


The  second  term  containing  r  and  the  first  partial  derivative  of  the 
concentration  present  no  problem  in  orthogonal  collocation  since  we 


have  the  discretization  for  all  partial  derivatives: 


N+2 


V6’ 


11+ 2 

j*1AijCA<rj't) 


(211) 


The  r  is  simply  replaced  by  the  at  each  collocation  point,  where 
the  x^  are  the  roots  of  the  approximating  polynomial  (i.e.  the  x^  or 
collocation  points  in  the  planar  cases) .  Ue  then  proceed  as  before 
to  solve  the  set  of  simultaneous  differential  equations.  We  reduce 
the  order  by  2  by  inserting  the  boundary  conditions. 


dCA 

ar 


N+l 

_  *  Bi,lcA(0,t)  4  Bi,N+2CAU#t)  +  HLBijCA(rjft) 


lAi,lcA(0,t)  4  Ai,N+2cA(1,t)  4  jB2AijcA(rj,t) 


Substituting  the  known  values  for  cA(0,t)  and  cA(l,t),  we  have 


dcA 

a r 


n+i 


11+1 


Bi,N+24j^2BijCA(rj't)  4  l“i,K+2'4*5”ijwA'* j 


+  l  AiAC%(xA,t) 

j»2 


(212) 


r 


(213) 


mp 

m 

|7> 

If 

m 

0 

f:> 
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At  this  point,  it  is  worth  considering  the  effect  of  choosing 
orthogonal  polynomials  other  than  the  Legendre  (y  *  6  «  0  in  equation 
(75))  polynomials.  Although  beyond  the  scope  of  this  chapter,  it  may 
be  shown  that  a  more  suitable  symmetry  for  spherical  diffusion  in¬ 
volving  the  differential  equations  used  here  is  one  in  which  v/e  choose 
•y  *  1  and  6  *  1/2. 

The  collocation  points  and  the  A. .  and  B. .  thus  are  different  from 

*  J  *»/ 

the  planar  diffusion  cases.  The  new  values  are  given  by  the  JCOBI 
and  DFOPR  subroutines  simply  by  changing  the  value  for  y  (AL)  and  6  (BE) 
in  the  main  program. 

Although  the  Legendre  polynomials  give  excellent  simulations  even 
in  this  case,  the  new  choice  of  y  and  6  is  even  better.  For  example, 
integration  of  the  equations  and  subsequent  calculation  of  the 
current  from 


I 


nr AD [A0] 

U-V 


nFADlA0)  N+l 

(L-It0)  lAl,N+2+^2AijCA(rjrt)1 


(214) 


leads  to  the  following  results  for  B*  ■  1  x  10  ^  in  equation  (213),  and 
time  point  t  *  0.03715: 


Iexact 

^simulated 

y  *  6  *  0  (Legendre)  0.937372 

y  -  1  ;  6  «  1/2  0.999958 


Six  internal  collocation  points  were  used  in  the  integration. 
The  additional  accuracy  can  become  important  in  dealing  with  very 
rapidly  changing  profiles. 

Spherical  symmetry  may  also  be  treated  by  recognizing  that  in 
such  a  case,  the  differential  equation  is  also  valid  for  all  radial 


directions,  i.e.  the  equation  is  valid  for  VC  as  well  as  It.  One  cay 
stake  the  variable  transformation 


u  *  r 


in  equation  (200). 
t;e  notice  that 


(215) 


dCA  _  d°A  du  r  =  ~  dCA 
dr  du  dr  du 


(216) 


I?  ’ 


+  4u 


(217) 


At  the  collocation  points,  equation  (211)  becomes 

fic*  dc. .  J„  d  CA  .  4u1/2  dcA 

Wt  „  2  du  du2  u1/  ~+g*  °u 


(218) 


with 


cA(0,t)  «  0 
cA(l,t)  *  cA(u,0)  *  1 


(219) 

(220) 


(221) 


Discretization  of  equation  (218)  with  insertion  of  the  boundary 
conditions  yields 


MW«  WT*  MT* 

Ht“  *  2  lAl,!H2'f;ji2AijCA(uj,t)1  +  4ui  lBi,N+2+^2BijCAluj't)] 

4ui1/2  N+l  6u£1/2 

+  lAi»n+2+jJ2AijcA(uj*t>I  "  u  W+e*  lAi,H+2 

N+l  N+l 

+j*2AijCA*Uj't>J  +  4ui  lBi,N+2+^2BijcA(uj#t)  J 


(222) 
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The  ik  are  simply  the  squares  of  the  values  of  the  collocation  points. 
This  set  of  equations  may  then  be  integrated  by  the  standard  nethods  to 
obtain  the  concentration  profiles.  This  set  of  equations  is  probably 
more  accurate  than  those  previously  described,  but  a  detailed  comparison 
has  not  been  made.  The  polynomials  used  here  again  should  be  those 
vhere  y  *  1  and  6  *  1/2. 

Example  5:  Chronopotentioir.etry  -  Planar  Electrodes 

(a)  Simple  Reversible  Charge  Transfer 
ne~ 

A  B  (223) 


The  diffusion  equations  and  boundary  conditions  describing  this 
experiment  are 


iz 


D  £-141 


(224) 


6  IB] 
6T 


B 


6HB] 

6X2 


(225) 


^X,0  *  W«,t  * 


(226) 


IB) 


X,0 


0 


(227) 


D  -  Da  .  Dfi 


(220 

(229! 

(230 

(231 


Insertion  of  the  usual  dimensionless  variables  for  concentration  and 
distance,  and  the  new  time  variable  t  *  T/t,  t  the  transition  time,  yields: 


IA#J6ca  dlA*U2cA 


x  fit 


L2fix2 


(232) 


IA*J6cb  DtA*]6  cb 


(233) 


T  fit 


L26x3 


and  after  simplifying  and  letting  8 


ffA  .  6  5  CA 


fit  fix 


(230 


*C3  _  £  6  CB 
fit  fix2 


(235) 


The  boundary  conditions  (226-228)  are  treated  similarly  and  become: 


cA(x,0)  *  1  cA(l,t)  *  1 

cB(x,0)  *  0  cB(l,t)  *  0 


K)  .  /M 

Vfix  /X1B0  {&* 


Discretization  of  equations  (234)  and  (235)  yields 


(236) 


(237) 


(238) 


dc.  N+2 

at-  *  b  1  B  .c  (x  ,t) 

dt  x  j-1  i#3  *  3 


(239) 


H+2 

8  I  B.  jC. (x . , t) 


(240) 


Expanding  (239)  and  (240)  partially,  we  have 
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6 1 Bi ,1CA (0 ' +Bi  ,H+2CA (1 ' +^2Bi ,  jCA (X  j ' t}  1 


(241) 


dt  x.  *  elBi#lCB(0#t)+Bi,N42CB(1,t>+^2Bi,jCB(xj#t)1 


(242) 


Insertion  of  the  boundary  conditions  (236)  and  (237)  yields 


6lBi,lCA(0#t)  +  Bi,N+2  +  .^xBi ^  jCA (x j " 5 


(243) 


MBiflcB(0,t)  4  Z*±lfy<*y t)3 


(244) 


We  need  an  expression  for  cA(0,t)  and  cB(0,t)  explicitly.  From  (92) 
applied  to  the  boundary  flux  condition  (238) ,  v;e  have 


x.-O  j 


N+2 

E  A-  .,c  (x.,t) 
j«=l  A*3  A  3 


(245) 


Xj*0 


M+2 

I  A,  .,C  (x,,t) 
A»3  ®  3 


(246) 


Now,  equation  (230)  is  also  trade  dimensionless  as  follows.  Since: 


w. 


1  vF/x.-o 


we  have 


>»P1A*I  (dCp\ 

L  Vdx  /Xl-( 


nrArA»jD1/2fi1/2 


Combining  this  expression  with  (245)  and  (246),  we  have 


(247) 


(248) 


cA(0,t)  «  Q  -  E  a.c.(x.,t) 

A  -U2  3  A  3 


(245) 


CB(0,t)  «  R  -  *fylXyt) 
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where 


0 


*1,1  [nFA[A*]DV*  Al'N+2j 

_J_  [ _ i.lL-,1 

*1,1  j_nFA  l A®  ]  D^B**  J  ,  and 


a  -  ^ 

2  *1,1  ' 

Substitution  of  (240)  and  (250)  back  into  (243)  and  (244),  we  have 


(251) 


(252) 


(253) 


dc. 


N+l 


•rrr  *  S.  +  6  '  r  b;  .c  (x.,t)  ,  and 

dt  i  j,2  i,3  A  3 


(254) 


dc 

It 


B 


N+l 


*=  T.  +  B  I  b.  .c  (x.,t)  ,  where 

A  j«2  1,3  B  3 

(255) 

BlBi,lQ  +  Bi,N+21  ' 

(256) 

®  Bi  1R  '  anfl 

(257) 

“  ~Bi,laj  +  Bi,j 

(258) 

Simultaneous  solution  of  the  2N  differential  equations  by  the  inte¬ 
gration  subroutines  described  previously  give  accurate  and  fast 
approximations  to  the  concentration  profiles  of  A  and  B  as  a  function 
of  tine. 

The  desired  response,  potential  as  a  function  of  time,  is  given 
then  by  equation  (231)  in  the  form 


RT 

nF 


(Cp ( 0 , t)\ 

c^Ttr)  < 


(259) 


with  cfi(0,t)  and  c^(0«t)  being  supplied  at  each  time  integration  step 
by  equations  (254)  and  (255). 


I 

| 

if 

m 

a;] 
1 
*  i 


% 
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(b)  EC?  „  Mechanism 

T6V 

The  extension  of  the  mathematics  to  include  kinetic  reactions  is 


(270) 


5  (il£L) 

C\6X  A 


Introducing  the  dimensionless  parameters  for  concentration,  time, 
and  distance,  the  following  simplified  equations  are  readily  obtained: 


(271) 


A (14X)  a(cb-Kcc) 


(272) 


6cp  6  cr 

—  .  e— r-  +  A  (l+K)  A  (cR-Xc_) 
fit  fix2  B  C 

cA(x,0)  «  cAU,t)  *  1 


(273) 


(274) 


cB(x,0)  «  cfi(l,t)  «  0 


(275) 


cc(x,0)  *  cc(l,t)  «  0 


(£)  --(£) 

\6x  Jx  \6x  yx  t 


f6cc>) 

W" 


(276) 


(277) 


(278) 


Equations  (271-273)  are  discretized  (93)  to  yield,  after  partial 
expansion  and  substitution  of  boundary  conditions  (274-276) , 


dcA(t) 


8(Bi,icA(°,t)  4  b1#n+2  4  ^^(Xj'tn 


(279) 


dcB(t)  N41 

-ar-  •  BrBi,icB(°,t)  4  - 


A(l4K)‘i(cB(x1,t)  •  Kcc(x1,t)J 


(280) 
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dcc(t) 


«  BIB,  .C  (0,t)  ♦  £  B.  ^C_ (x. rt)l  + 

Xi  1,1  c  j*2  i#3  C  3 


X(l+K)'itcB(xi,t)  -  XccUift)l 


(281) 


The  flux  relation  (277)  is  identical  to  the  simple  reversible  charge 
transfer  case,  and  is  discretized  in  precisely  the  same  nanner  to 
yield  equations  (249)  and  (250) .  The  flux  relation  for  species  C  is 
equation  (278),  and  is  discretized  as  follows: 


11+2 

*  0  *  £  A,  .c_(x.,t) 


j*l  C  D 

N+l 

Ai,icct0,t,+Ai,:i+2cc(1't)+^2Ai,  jcc(xj#t) 


A.  -Cr<(0,t)  +  £  A.  .c  (s:  ,t) 

1#1  C  j«2  -1'3  C  3 


(282) 


cr(0,t)  «  £  a,c_(x. ,t) 
c  c  J 


where  a^  was  defined  by  (253) . 


Fufcstitution  of  (249),  (250)  and  (283)  into  (279),  (280),  and  (281) 
gives 


(283) 


dCA 

2RT  * 

S,  +  $ 

;{i 

dCB  . 

ar  ■ 

Ti  +ji 

dcr 

11+1 

3T  * 

Xj 

8  E  b. 
j-2  x 

M+l 
£  l 
j*2 


(284) 


(285)  Wt 

•if 

* 

■i 

(28€)  | 


As  before,  simultaneous  solution  of  equations  (284)  -  (286)  provides 
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the  tine  dependent  concentration  profiles.  The  chronopotentionetric 
response  for  this  mechanism  is  still  given  by  equation  (259)  with  the 
term  in  parentheses  being  supplied  by  equations  (249)  end  (250) .  The 
concentration  terms  in  (249)  and  (250)  however,  are  now  calculated 
numerically  from  equations  (284)  -  (286). 

The  elements  of  the  Jacobian  matrix  are  simply  the  derivatives  of 


equations  (254-255)  or  (284-286)  with  respect  to  each  x^.  The  matrix 


elements  of  the  Jacobian  for  the  simple  reversible  electron  transfer 
mechanism  then  are  given  by: 


i.i  SKj  \  at  j 


for  species  A 


x. 

i 


and 


'i,j  sxj  \  at  J 


for  species  B 


These  are  given  explicitly  by: 


S  (  N+1 

ri,j  “  ^Si  +  £2hi,ic*ixy 


t) 


for  species  A,  and 
J 


*  (  N+1  \ 

i.i  *  Sj(Ti  +  j^bi.iCB<Xi't,J 


for  species  B,  or  finally, 

’i.i  *  bi.i 


1^1 

jj  4  *  b,  j  *  — +  b 


i#i 


for  both  species. 


(287) 


(288) 


(289) 


(290) 


(291) 


0 

V  'i 


m 

so 

•yf 
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Similarily ,  the  matrix  elements  of  the  Jacobian  for  the  EC 
mechanism  are  given  by: 


J.  .  *  equation  (291) 


for  species  A,  and 

h.i  -  ~B^,j  4  ‘  «(i«)-1cB(*1,t)slj 

for  species  B,  where  6 _  is  the  Kronecker  delta,  and 

Jj  «s  *  ■  +  B.  .  -  XKd+Kj'^Cx.  ,t)6. 


rev 


(292) 


(293) 


'i,j 


*1,1 


'C i uij 


(294 


for  species  C. 

Listed  below  are  the  discretized  equations  for  several  other 
commonly  occurring  electrochemical  mechanisms.  As  above,  it  is 
assumed  that  the  diffusion  is  to  a  planar  electrode  in  quiet  solution, 
and  species  A  is  the  only  electroactive  species  present  initially. 

(c)  Catalytic  Mechanism 
ne~  k 

A  B  — >  A  +  C 


dcA 

a r 


dcB 

ar 


N+l 

Si  +  M  bi  .Cft(X.,t) 
1  j«2  1,3  8  3 


N+l 

Ti  + 


ktcBCXj,t) 


ktCgUjit) 


(295) 


(296) 


(297) 


For  the  E-t  profile,  equation  (231)  is  used  along  with  equations 
(249)  and  (250) ,  the  unknown  concentration  terns  being  furnished  by 
the  simultaneous  solution  of  the  2N  equations  (296)  and  (297). 

(d)  Dimerization 
ne" 

a  b 


(298 


/m 

S-VjpS; 

X 


4 


Si 
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2B 


k 


1 


C 


(229) 


Me  let 


X  «  ( IA*]K1  +  k^1)T,  and 
^-1 

K  *  WWl 


Then 


Si 


N+l 
e  i  b 
j«2 


ffE 

dt 


T. 


1 


+ 


N+l 

Mb.  .c_(x.,t)-X(l+K) 
j*2  1,3  B  3 


1(cB(xi,t)2-Kc£;(xi,t)) 


N+l 

I  b,  jC_ (x . , t)  +  X (1+K) 
j*2  1,3  C  3 


1(cB(xi,t)2-Kcc(xi,t)) 


(300) 

(301) 

(302) 

(303) 

(304) 


Once  more,  equations  (231),  (249),  and  (250)  are  used  for  the  E-t 
profiles.  Note  that  in  this  case,  3N  equations  must  be  solved 
because  c£(0,t)  depends  on  the  concentration  of  the  C  species. 


(e)  Second  Order  Reaction 


X  and  X  are  defined  as  in  equations  (300)  and  (301) . 


dcA 

IT 


N+l 

*s .  +  e  i  b 

j-2 


i^cA(Xj,t)-X(l+K)~*(cA(Xi,t)cB(xi»t)- 


(305) 


Kcc(x.,t)) 


(306) 
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dcR  N+l 

an.  ‘  Ti +  6^2bi.jcB(xj-t> +  X(1+K1  ‘cAtei't»e#<*i't> 


KoC  <2:i '  t) ) 


(307) 


Ssl  . 

dt  K 

Ki 


Bib.  .,cr(x.,t)  +  Hl+K)‘1(c.(x.(t)cR|s.  ,t)-Xcr(x.  ,t))  (308) 

j*2  1,3  C  3  A  i  B  i  Ci 


F.  Cyclic  Voltammetry  -  Planar  Diffusion 


nieker  and  Speiser  [7]  have  derived  the  discretized  equations 
necessary  to  simulate  cyclic  voltannetric  responses  by  orthogonal 
collocation,  Tney  are  presented  in  this  section.  The  only  dif¬ 
ference  in  the  mathematical  treatment  for  this  case  and  chronoamp- 
eronetry  is  the  statement  of  the  surface  concentrations  of  each  species, 
which  of  course  must  be  modified  to  reflect  the  nev;  potential  program.  It 
should  be  noted  that  the  computer  time  expended  in  these  cases  is  approx¬ 
imately  the  same  as  Feldberg's  finite  difference  approach. 

(a)  Discretization  of  Cyclic  Voltanmetric  Experiments  (7) 

ne” 

1.  Simple  electron  transfer  A  B 


From,  the  original  boundary  conditions  and  differential  equations, 
we  have  the  following  dimensionless  equations: 


«cR  62cR 

—  -  B - 1 

«t  6x* 


cA(x,0)  «  1,  cB(x,0)  »  0 


cA(l,t)  *  1,  cB(l,t)  •  0 
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*CA  6CB 

*="  *.o  *  '  ="  x-0 


(313) 


ca(0,t) 

c^7o7t)  *  ®A/B  SX(t) 


(314) 


with  ®A/B  “  e*PlnF/RT  <ea/B  “  E8tartn  t315) 

S^(t)  *  |exp(-nFvt/RT)  «  exp (-at)  t  S  t^  (316) 

|exp(at-2atx)  tT  5  t  <  2tA  (317) 

*  potential  sweep  switching  tine 

a  «  nFv/RT  (313) 

dE 

v  **  cTt'  the  8can  rate  (315) 

EA/B  *  the  formal  potential  of  the  A/B  couple 

Estart  *  the  8tart*n9  potential  of  the  scan 

8  -  D/aL2  (320) 


The  differential  equations  (308)  and  (310)  are  discretized  by 
the  approximation  (93)  at  the  collocation  points  x^. 

N+2 

-BIB.  .c.  (x . ,  t)  (321) 

xA  j-1  1,3  A  3 


dcB 

ar 

Expanding  out  the  two  boundary  terms  at  x  -  0  and  x  *  1,  i.e.  xi  -  xQ 
end  x^  -  Xjj+2»  w®  have 


H+2 

-BIB 


i,jCB(xjrt) 


(322) 


6lBi#lCA(0,t)+Bi/N+2cA<^t)  ♦  IB i#jCA(Xj,t) 


(323) 
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at  *  Bl£i,lCB(0,t)+Bi,IJ+2CD(1,t)+.Z^Di, jcB(xj#t) 

X, 


(324) 


£?ut  stituting  the  boundary  conditions  .3  and  .4  yield 


ot  L 


6 lEi , 1CA (0 ' +Bi ,N+2+jZ 2Bi , jCA (x  j ' t} 3 


(325) 


ZtT  *  BlBi#lCB<0,t)+ABirjCB(xj'1 


(326) 


An  explicit  value  for  c%(0,t)  and  c,-(0,t)  rav  be  derived  exactly  as  v/as 

r.  x> 

done  for  previous  equations  by  usinq  the  boundary  condition  (313) .  Ihe 
interr..ef'  iate  result 

"1 , 1CA ( 0 ' 15  +A1  ,N+2+  Z  A1 , jCA j ' t}  *  *IAl,lcB(0,t> 


iH  1 

I  A.  •  c_  (x . ,  t)  J 

a  J 


(327) 


is  nov?  modified  to  include  the  tine  dependence  of  potential  and  to 
eliminate  one  of  the  unknowns  cA(0,t)  or  cB(0,t)  by  the  insertion  of 
the  boundarv  condition  (314). 


The  result  is 


cA(0,t) 


eA/BSX (t) 


A,  Z  A.  4lc.  (x.,t)+cn(x.,t)]  (328 


Aifi(l+eA/BSx(t)]  “1 ,11+2  1 ,  j 1  A  '  j  B  j 1 


CB(0,t)  • 


"l,llA  °A/B  a 


Al,M+2+^2Al,jlcA(xj,t)','CB(Xj't)1  (329 


Insertion  of  these  explicit  values  of  cA(0,t)  and  cB(0,t)  into  (325) 
and  (320  yield  the  final  set  of  N  differential  equations  to  be  solved 
simultaneously  for  the  N  values  of  cA(x^,t)  and  c&(Xyt)  (at  the 
collocation  poir.cs)  * 
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dc; 

It 


B 


,e*,rA  (t) 


•l!lll+eA/BSA  (t)]  L  1  'N+2+  *  A1 ,  j  1 CA (X  j ' t}  +CB  (X  j  ' t}  ]  J 


N+l 


(330) 


dc 

it 


B 


Xi 


B.  x 

Al,lll+6A/BrX(t) 


r  N+l  *1 

T  [A1  ,N+2+  . ^2A1 ,  j  1 CA  (X  j  ' 1 5  +CB (X  j  ' t} 


N+l 

+  IB,  .C-.  (3C  .  ,  t) 

jc2  x*3  B  j 


(331) 


The  concentration  profiles  are  fully  developed  fror.  the  solution 
of  these  equations. 

The  current  is  given  by 


6c, 


i  *  nFA'D 


dc. 


,n  -nFA'/oS  (A«3sr 


x^O 


(332) 


where  A'  is  the  electrode  area. 


dc. 


The  discretization  of 


ax 


follows  fror. 


x^O 


dc 

it 


A 


Xj*0 


N+2 

j«l  A»J  A  3 


(333) 


Expansion  of  (333) ,  insertion  of  the  boundary  conditions  (312)  and  (314) 

with  insertion  of  the  result  back  into  (332)  yields 

nFA'/DSa  [A° ]  r  N+l  N+l  T 

1  *  1+eA/BSX(t)  [Al»N+2+j^2AlrjCA(Xj't)“6A/BSX(t)  ^l,  jCB(Xj 't}J  (33< 


Since  we  know  the  c.(x.,#t)  and  C-(x.,t)  fron  the  solution  of  (330)  and 

A  D  B  j 

(331) ,  the  current  as  a  function  of  tire  is  irmediately  found.  Since 
the  real  titae  T  *  at  is  related  to  the  potential  by 


*  -  Wt  +  **  **  *  *  *  0  0  it  £7X 


(335) 


I 
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E  *  EStart  +  v(2VTI 

Tx  5  T  <  2TX 

(336) 

or  dimensionless  tine  by 

E  •  E£tart  +  vat  5 

'  *tSt\ 

(337) 

E  "  Estart  +  va(2Vt) 

t^  ^  t  ^  2tx 

(33S) 

The  spectroelectrochenical 

absorbance  response  A  is  given  for  each 

species  A  and  B  for  cyclic 

voltammetry  by  the  quadratures 

N+2 

A,  «  y  £  V, c 
'•  i-i  1  Ai 

(339) 

N+2 

ab  *  y  I  v. c_ 

B  i=l  1  Bi 

(349) 

where  the  are  given  by  the  weight  vector  formula  solved.  In  the 
programs,  the  A  are  given  at  any  tine  t  simply  by  using  the  c-(x.,t) 

n  1 

or  cD(xi,t)  in  (339)  or  (340)  respectively,  that  correspond  to  that 
tine  t. 


2.  Reversible  electron  transfer  followed  by  irreversible  first 
order  chemical  reaction: 


-ne  k 
+ne” 


C  <ECI!tt> 


The  dimensionless  partial  differential  equation  for  species  B  is  now 

«cf 


fit 


-B  ‘S.  ae 


6x‘ 


B 


(341) 


v’here 


© 

inr 

i 


8 


a  »  k/a 


(342) 
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9y  inspection  of  the  equations  for  the  reversible  charge  transfer 
redian  it  is  seen  that  only  the  equations  describing  the  B  species 
need  to  be  changed  to  fit  the  ECIRR  case*  The  changes  yield  these  new 
equations  for  B: 


dc 

St" 


B 


N+J 


BlBi,lCB(0't)+Bi,N+2CBat}+jf2Bi,jCB(xj»t)}“0fCB(Xj»t)  (343) 


dCj 

dt 


xi 


Al,lI1*®A/BSX(t)]  [Al»N+2%LAl^ICA(Xj,t>+CB(: 


N+l 
■  l  I 
j«2 


Cj't)3] 


N+l 


+  IB.  jC—  (X  .  r  t) 


jt=2 


i  #  j  B  j 


-  ocB(xi#t) 


(344) 


The  current  is  exactly  the  same  in  the  forward  direction  as  in  the  pre¬ 
vious  case  since  the  kinetic  tern  does  not  affect  the  flux  of  species 
A  to  the  electrode.  The  spectroelectrochenical  absorbance  is  also 
identical  for  species  A#  but  is  changed  for  species  E. 


3.  A 


“V 


B 


+nle 


B 


H+  +  C 


H+  +  A 


HA 


-n2e 


HA  v  D 


+n«>e 


k^  (first  or  pseudo  first  order) 
k2  (second  order) 

EHA+/D 


(345) 


D  — ♦  2H+  +  C 


k3  (first  or  pseudo  first  order) 


The  results  are  given  below  for  this  complicated  mechanism.  The 
reader  is  referred  to  the  original  literature  [7]  for  the  details  of 
the  derivation. 
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fS 

dt 


dc 

dF 


B 


i  b.  r  tj+i  1 

Ki  -  4  [^v^,CA<Xj’t,+CB  ^’’J 

N+l  ) 

+  Bi .N+2* ^ jBi , jCA <X j ' 1  *  j  '  a2CA(Xi't)V(Xi't)  < 

Xi  *  4  Aa.iU^rTtTT 

N+l  |  , 

+  41Bi.jcB(x:'t,j  *°icB<xi't>  ‘ 


dcc 

ar 


8  N+l 


N+l 


Xi 


g[-  -i Ai  £  A.  .c  (x.,t)+  Z  B.  .c  (x.,t)] 
Al,l  j«=2  X'j  C  3  j-2  A'3  C  3 


+  a1cB(xi,t)+a3cD(xi,t) 


(343) 


^CH+] 

<St 


B.  .  N+l  N+l 

-  B[-  -ixi  I  A.  .C„+(X.,t)+  Z  B.  .c  +(x.,t)] 
Dt.  *1,1  j-2  1>3  H  3  j-2  1,3  “  3 

+  a1cfi  (xA ,  t) +2a3cD  (Xi ,  t)  -a2cA  (x^, ,  t )  cfl+  (xA ,  t) 


(342) 


dCHA+ 

dt 


A  Bi,ieiIAVDSl(t) 

j  Alflll+eHA+/DS^lt  ^ 


N+l 

j^2Al,j[cHA+(xj 


'tl+CglXjft)] 


N+l 

+  j^2Bi#jCHA+(Xj,t) 


+  02CA(Xirt)CK+(Xi#t) 


(350) 


B1 . 1 


N+l 

1LAi.j|CHA+<Xj't,+CI>(!tj't)! 

J 


N+l  < 

+  ji2Bi>3CD<X3't>l  '  “S'D111!'11 


(351) 


In  the  above  equations, 
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J 


°1  *  kl/a  (352) 

a2  *  k2tA*1/a  (353) 


°3 


(354) 


6HA+/D  *  exp  |r,r/RT(EAAVD  -  ^tart'1 


(355) 


The  current  is  given  by 

dc.  dc„,i 

i  *  FA  D6a  [A°]  [n.  +n,  -3-  J 

1  8*  Xj-9  2  ox  Xj-0 


(356) 


or 


{  n.  11+ 1 

FA  D6alA#jji+e^sS^t)  [A1>N+2+^I^A1^CA(xjft)- 


A/S  A 
N+l 

~  ea/np%(t>  *  A,  .C-,(X.,t))  + 


j-2 
N+i 

“6HA+/DSA{t)  3f2Al#jCD(xj#t)1| 


H+l 

I  £  A,  4c«»+(x^»t) 


T  i+W/DSA(t)  lj«2  I1A+  " 


(357) 


(b)  Optimization  of  the  Dimensionless  Parameter  $ 


The  dimensionless  parameter  6  that  is  introduced  into  the  col¬ 
location  equations  must  be  chosen  such  that  (a)  it  is  not  so  lov;  as 
tp  distort  the  simulations,  and  (b)  it  is  not  so  high  as  to  cause 
oscillation  at  a  collocation  point.  Since  the  parameter  $  is  a  function 
of  the  dimensionless  rate  constant,  simply  choosing  a  fixed  value  will 
not  give  the  best  results  for  all  possible  values  of  rate  constants. 
Speiser  (18]  has  recognized  this  fact  and  has  developed  an  analysis 
of  the  problem  based  on  stability  criterion  for  the  particular  integration 


method  used.  An  optimal  value  of  6  then  is  calculated  in  the  program 
itself. 

A  numerical  integration  is  "stable"  if  the  difference  between  the 
true  and  approximate  solution  decreases ,  while  the  condition 


6f (x,y) 
6dy 


<  0 


(353) 


where 

f  (x,y) 
holds . 


dv 

dx 


(355) 


If  v;e  are  referring  to  a  system  of  II  ordinary  differential  equations, 
then  y,  f(x,y),  and  - represent  column  vectors  with  N  components. 
For  example,  in  the  EC___  mechanism,  the  components  of  the  vector 


^dc(xi,t)^ 

6c(x,t) 


(360) 


are  the  partial  differentials  of  equation  (330)  and  (344)  and  are 


IcTxTTty 


Bi,leA/BSl(t) 

Al,lli;eA/BElf0j 


Al,i+Bi,i 


2,  ...N+l  (361) 


B 


if  1 


+  e  i-2,  ...N+l  (362) 


It  is  noted  that  for  the  values  of  the  A.  .  and  B.  4  found  in 

+ij 

program  DFOPR  (see  Appendix),  the  expressions  (362)  and  (363)  are 
always  negative,  and  hence  always  stable. 

For  the  Hammings  predictor-corrector  method,  such  as  that  found 


in  IBM's  SSP  library,  Speiser  points  out  that  a  particular  stability 
criterion  for  our  problems  is 


and  are  found  to  be: 
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0.65 


I  El.ieHf.VDS1)'^t) 


Al,l(l+6HAV3SX(t,J  Al'i+!>t-i 

0.65-ho2cA (x^t) 


Ai  <+B<  < 

^2  **i  l*i 


#  i*2|  . • .N+l 


,  i«2,  ...N+l  (366) 


Thus  using  the  inequalities  defined  in  (365)  and  (365) ,  and  in 

(366)  one  nay  optimize  6  such  that  the  integrations  approach  the  highest 

possible  accuracy  without  going  into  oscillation.  Speiser  notes  that 

negative  values  of  6  cannot  be  used  since  the  current  is  proportional 
1/2 

to  B  ,  so  negative  values  of  6  must  be  eliminated  for  a  given  set 
of  kinetic  parameters  by  decreasing  h  until  0  becomes  positive. 
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APPENDIX 


Subroutine  JCOBI:  Calculates  roots  and  derivatives  of  a  polynomial . 

*Usage:  CALL  JCOBI  (ND,  N,  NO,  Nl,  ALFA,  BETA,  DIFl,  DIP2,  DIF3,  ROOT) 

Input  parameters: 

Integer  ND:  The  dimension  of  vectors  Diri,  DIF2,  DIF3,  and  ROOT. 

These  vectors  should  be  dimensioned  in  a  DIMENSION 
statement  in  the  calling  program.  ND  should  be  at 
least  as  large  as  the  total  number  of  interpolation 
points. 

Integer  N:  The  degree  of  the  Jacobi  polynomial  (i.e.  the  number  of 
interior  interpolation  points) . 

Integer  NO:  Control  index.  If  Nl*l,  x*0  is  used  as  a  collocation 
point;  if  N0=0,  x«=0  is  not  used. 

Integer  Nl:  Control  index.  If  Nl*l,  x«l  is  used  as  a  collocation 
point;  if  N1«0,  x«l  is  not  used.  Normally,  boti.  NO  an 
Nl  are  set  to  1. 

Real  ALFA, 

BETA:  Determines  the  symmetry  of  polynomial  chosen.  (See 

equation  (75  ) ) .  The  Legendre  polynomials  described 
herein  have  ALFA*BETA»0  for  linear  diffusion  problems. 

Output  parameters: 

Real  Array 

Root:  One  dimensional  vector  containing  on  exit  of  the  sub¬ 

routine  the  N+N0+N1  roots  or  zeros  of  the  chosen 
polynomial . 

Real  Array 

DIFl,  DIF2, 

and  DIF3:  Three  one  dimensional  arrays  containing  on  exit  from 

JCOBI  the  first,  second,  and  third  derivatives  of  the 
polynomials  at  the  roots  of  the  polynomial. 


* 

Please  note  that  in  order  to  retain  Villaden's  original  program,  ALFA 
and  BETA  are  used  here  for  y  and  6  in  the  text  (equation  (75) ) . 


1 

SUBROUTINE  UCOBKND.N.NO.NI ,  AL.8E .OIF  1 .Dlf2.01F3.R00T) 

3 

IMPLICIT  REAL *8(A-H, 0*2 ) 

3 

DIMENSION  OIF  UNO) ,0IF2(ND) ,DIF3(ND) .ROOT(NO) 

4 

AB*AL+BE 

8 

AD-BE-AL 

S 

AP*8E  *AL 

7 

OIF 1( 1 )•( A0/( AB*2  >♦ 1 3/2 

• 

D1F2( 1 )*0. 

• 

IF(N.LT .2)60  TO  IS 

10 

00  10  1*2. N 

11 

21*1-1 

12 

2*AB*2*21 

13 

OIF  11 1 )*(AB*A0/2/(2*2)*1 )/2 

14 

1F( I . NE .2)00  TO  11 

IS 

01F2(I)*(AB*AP*21)/2/2/(2*1) 

IS 

CO  TO  10 

17 

11 

2*2*2 

18 

V*21*(AB*21) 

19 

V*V* ( AP*V  ) 

20 

D I F  2 ( I )*V/2/(Z- 1 ) 

21 

10 

CONTINUE 

22 

15 

X*0. 

23 

00  20  1*1, N 

24 

25 

X0*0. 

25 

XN*1 . 

26 

X01-0. 

27 

XN1*0. 

28 

DO  30  U-1.N 

29 

XP*(DIF1(U)-X)*XN-01F2(J)*XD 

30 

XP1»(DIF1<d)-X)*XN1-0IF2(J>*XD1-XN 

31 

N1*1 

32 

X0*XN 

33 

X01-XN1 

34 

XN-XP 

35 

30 

XN1-XP1 

36 

2C*1 . 

37 

2-XN/XN1 

38 

IF( I . EO. 1 )G0  TO  21 

39 

DO  22  U-2.I 

40 

22 

2C*2C-2/( X-ROOT ( 0- 1 ) ) 

41 

21 

Z*Z/ZC 

42 

X*X-2 

43 

I F(0A8S( 2 ).CT. 1.0-09)60  TO  25 

44 

ROOT( I  )»X 

45 

X*X*.0001 

46 

20 

CONTINUE 

4’ 

NT*N*N0*N1 

48 

IF(NO. 60.0)00  TO  35 

49 

00  31  1*1, N 

SO 

U*N*  1-1 

SI 

31 

R00T(J41)*R00T(J) 

S2 

ROOT( 1 )*0 

S3 

35 

1F(N1 .20. 1 IROOT (NT )• 1 , 

84 

00  40  1*1, NT 

86 

X*R00T( 1 ) 

86 

OIF  1(1 1*1. 

87 

01F2( 1 1*0. 

88 

01 F3( 1 )*0. 

89 

00  40  U*1,NT 

90 

IF(0.C0.1)00  TO  40 

81 

V*X-ROOT(d) 

82 

01F3( l )*V*DIF3( I )*  9*01 F  2(1  ) 

83 

0IF2O  >*Y*01F2( 1 )+2*DIF  1( 1 ) 

84 

OIF  1(1 )*V*01F1(I ) 

9S 

40 

CONTINUE 

86 

RETURN 

87 

CNO 

Subroutine  DFOPP.:  Calculates  the  A^,  ,  and  Qi  interpolation  and 

quadrature  weights. 

Usages  CALL  DFOPR  <ND,  N,  NO,  Ml,  I  ID,  DIFl,  DIF2,  D1F3,  ROOT) 


Input  parameters: 


Integer  ND, 

N,  NO,  Nl:  Sane  as  in  JCOBI 

Integer  I :  Index  of  the  root  at  which  A .  . ,  B . . ,  or  Q,  is  being 

evaluated.  i3  13  1 


Integer  ID:  Control  index.  ID=1  causes  calculation  of  A..,;  ID=2 
is  for  the  B^;  ID«3  is  for  the  Q^.  13 

Real  Array 
ROOT,  DIFl, 

DIF2,  DIF3:  The  vectors  calculated  in  JCOBI. 


Output  parameters: 

Real  Array 

VEC:  The  vector  containing  the  A..,  B,,»  or  Q.  at  the  given 

root  upon  exit.  13  13  A 


1 

SUBROUTINE  OFOPR(ND.N.NO.N1,1.ID.DIF1.D1F2.DIF3.ROOT.VECT) 

2 

IMPLICIT  REAL*8(A-H.0-Z) 

3 

DIMENSION  DIF1(ND),01F2(ND),D!F3(N0) .ROOT (NO) .VECT(ND) 

4 

NT«N*N0*N1 

S 

IFUO.EO  3)00  70  10 

c 

00  20  0*1, NT 

7 

IF (U .NE . I  )G0  TO  21 

• 

1F( IO.NE . 1 )C0  TO  5 

9 

VECT( I )>0IF2( I )/DlF 1(1 )/2 

10 

CO  TO  20 

11 

5 

VECT(I)-01F3(I)/0IF1(l)/3 

12 

CO  TO  20 

13 

21 

V-ROOT(I)-ROOT(d) 

14 

VECT (J)»01F1(I)/DIF1(J)/V 

15 

IF( 10  E0.2)VECT(d)-VECT(J)«(01F2(I)/DIF1(I)-2/Y) 

IS 

20 

CONTINUE 

17 

CO  TO  50 

19 

10 

V«0. 

19 

DO  25  J>1,NT 

20 

X-ROOT(J) 

21 

ax-xmi-x) 

22 

IF(N0  E0.0)AX«AX/X/X 

23 

1F(N1.I0.0)AX»AX/<1*X)/(1-X) 

24 

VECT(d)»AX/OIF 1(d)**2 

25 

25 

V»V*VECT(d) 

2fi 

00  CO  J-1.NT 

27 

CO 

VECT(J)-VECT(d)/Y 

28 

50 

RETURN 

29 

END 

Subroutine  STIFF3:  Perfoms  the  integration  of  single  or  coupled  systems 

of  stiff  differential  equations  by  the  ISI3  method. 

Osage:  CALL  STIFF 3  (N,  ND,  HPRINT,  FUN,  DFUH,  OUT,  XO,  XI,  HO,  EPS, 

V,  y,  YOLD,  YOLD1,  IP,  YA,  YKl,  YK2,  YK3,  DF,  DFOLD,  F,  FOLD) 

IP  is  an  internally  used  integer  vector  and  should  be  dimensioned  as 
such  in  the  calling  program  with  dimension  ND.  The  following  are  real 
rectors  of  dimension  1ID  and  should  also  be  dimensioned  in  the  calling 
progran  as  such:  If,  Y,  YOLD,  YOLD1,  YA,  YKl,  YK2,  YX3,  F,  FOLD.  All 
but  V  and  Y  are  internally  used  only.  DF  and  DFOLD  are  internally  used 
arrays  of  dimension  (HDXND)  and  should  be  designated  as  such  also. 

FUN,  DFUN,  and  OUT  are  external  subroutines  called  from  STIFF3  and  should 
appear  in  an  EXTERNAL  statement  in  the  calling  program. 


Input  parameters: 


Integer  N:  Number  of  equations  to  be  integrated. 

Integer  ND:  As  in  JCOBI  -  main  progran  array  dimension. 

Integer 

HPRINT:  Printing  interval.  The  solution  is  printed  at  every 

NPRINT  step.  The  solution  is  always  printed  at  Xl. 


Real  XO, 

Xl:  The  limits  of  tine  over  which  the  differential  equa¬ 

tions  are  integrated. 

Real  KO:  Suggested  initial  half-steplength  for  integration.  On 

exit,  HO  contains  suggested  value  for  half-steplength 
for  continued  integration  beyond  Xl. 


Real  Vector 

W  and  real 

EPS:  These  parameters  dete — tne  subsequent  steplengths  for 

integrations.  They  >•  part  of  the  automatic  step  size 
selection  part  of  the  program.  Suggested  values  are 
given  in  the  examples.  Further  information  is  given  in 
the  original  literature. 

Real  Vector 

Y;  Vector  of  concentration  solutions  at  the  collocation 

points.  Initially,  these  are  specified  as  the  tof  mdary 
conditions  (t*0)  for  each  species. 

6TIFF3  also  uses  three  internal  subroutines,  SIRK3,  LU  and  BACK.  The 
listings  for  theue  are  included. 


1 

2 

3 

4 


SUBROUTINE  SIRK3(N.NO .FUN. IPI V.F . V . VK1 . YK2. VK3.0F ,M) 

IMPLICIT  REAL»8<A-H.O-2) 

DIMENSION  F(NO),Y(ND),YK1(ND), YK2CND  J , ¥K3(ND) , IPI V(ND) , OF (NO. NO ) 
DIMENSION  R(4) 


5 

6 
7 
6 
9 

10 

11 

12 

13 

14 

15 

16 
17 
19 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 


OATA  A. R/. 43586652 15084S89D0. 1 .037609496131859D0,  .834930483852637 
1700. - .630202088724452300. - . 2423378912600452/ 

DO  5  I»1. N 
DO  6  J* 1 ,N 
DF(I.U)— M*A»DF(I.d) 

IF (DABS (DF (I,J)).LT.1,D-12)DF(I.J)«0. 

6  CONTINUE 

5  DF( I , I  )“DF (1.1  )♦ 1 . 

CALI  LU(ND.N.IPIV.DF) 

CALL  9ACK(N0,N, IPIV.DF.F) 

DO  8  1*1. N 
VK1 ( I )«H*F( I ) 

8  VK2(U«V(I>*.7500*YK«'I) 

CALL  FUN( YK2.F ) 

CALL  «ACK(ND.N, IPIV.DF.F) 

DO  9  I ■ 1 ,N 
YK2( I )*H*F( 1 ) 

V(I)»Y(I)*R(1)*YK1(I )*R( 2 )*YK2( 1 ) 

9  VK2(I)»R(3)»YK1(I)*R(4»*YK2<n 
CALL  BACK(NO,N. IPIV.DF . YK2 ) 

DO  10  I* 1 ,N 

10  V( I )•¥( I )*YK2( I ) 

RETURN 

ENO 


10 

11 

12 


2 


SUBROUTINE  ST1FF3(N.N0,NPR1NT .FUN.DFUN.0UT.X0.X1.H0.  CPS. R.Y. YOLO. 
1 YOL0 1 . 1 P . V A . Y« 1 . YK2 . YK3 . OF . DFOLD . F . FOLD ) 

IMPLICIT  REAL*9( A-H.O*Z ) 

DIMENSION  IP(NO),V(NO), VOLD( NO ) , V0LD1 (ND),YA(ND),VK1(ND),VK2(ND) 

DIMENSION  YK3(ND ) , W(ND ) , F (NO ) ,FOLD(ND ) .OF (NO ,N0 ) , OF OLD (NO. NO) 

IC0N*0 

NOUT*0 

X-XO 

H*HO 

IF  (XO  ♦  2.»H  .LT.XOCO  TO  1 

H*(X1*X)/2 

ICON*  1 


13 

14 

15 

16 
17 
IS 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 
.  97 

SB 

39 

40 

41 

42 

43 

44 

45 

46 

47 
49 
49 
80 
81 
82 

83 

84 
88 
86 
87 
89 
89 
60 
61 
62 

63 

64 

65 

66 
67 


1  IF  (ICON  .CO.  0  .AND .  X«4*H  ,CT.  X 1 )H«( X 1-X )/4 

CALL  FUN(V.F) 

CALL  OFUN(V.OF) 

IHA.-1 
DO  30  I • 1 ,N 
V0LD( I )*Y( I ) 

FOLD( I )*F ( I ) 

DO  30  d* 1 ,N 

30  DFOLD( I , d)*DF (1,0) 

37  CALL  SIRK3(N,ND,FUN, IP.f ,Y. YK1 .VK2.YK3.DF ,2*H) 

DO  35  I* 1 ,N 
VA( I  )*V( I  ) 

V { I )*V0L0(1 ) 

F ( I  )«FOLD( I ) 

DO  35  d*1.N 

35  DF( I , 0 )*DFOLD( I  ,d) 

39  IH4* IHA+ 1 

CALL  SIRK3(N.N0.FUN.IP,F,Y.YK1,YK2.VK3.DF.H) 

CALL  FUN(Y.F) 

CALL  DFUN(Y.DF) 

DO  40  1-1. N 

40  YOLDI(I)-YU) 

CALL  SIRK3(N.N0,FUN.IP.F.Y.VK1,YK2.YK3.DF,H) 

C*0. 

DO  41  1-1. N 

CS*W(I)»0ABS(VA(I)-Y(I))/(1.*DABS(Y(I))) 

IF(ES  GT.  C)C*ES 

41  CONTINUE 
O-E/EPS 

QA«(4.-Q)**.25 
I F ( 0  .LE.  1.)  GO  TO  49 
DO  48  I-1.N 
VA(I)-YOLOKI) 

F(l)-FOLDd) 

V(n*VOLD(I) 

DO  45  0*1. N 

45  OF (I ,U)*DFOLD(I .0) 

H*H/2 
I CON-O 
00  TO  36 

49  00  49  I-1.N 

48  V( I )-Y(I )♦( Y(l )-VA( I ) )/7 .00 

X*X*2*H 

0A- 1 . /(0A« 1,0*10) 

IF(OA  ,0T.  3. )0A*3. 

M-OA*H 
NOUT -NOUT ♦ 1 

IF ( (NOUT /NPB1NT )*NPBINT  .10.  NOUT  .OR.  ICON  .CO.  DCAU  OUT(X.Y  .1 
1HA.0A) 

IF ( ICON  .CO.  1)  00  TO  187 

H0‘ 

IF(»«2.*H  .IT.  XI)  00  TO  1 
00  TO  2 
167  RETURN 
6  NO 


1 

2 

3 

4 

5 
C 
7 
• 
9 

10 

11 

12 

13 

14 
19 
1C 

17 

18 

19 

20 
21 
22 

23 

24 

25 
2C 

27 

28 

29 

30 


SUBROUTINE  LUIND.N. IPIV.A) 
IMPLICIT  REAL *8 ( A -H, 0*2 ) 
DIMENSION  IPIV(ND).A(NO.ND) 
IP1V(N)«N 
N1«N-1 

00  10  1»1.N1 
X«A(1,I) 

1F(X.LT.0.)X«-X 

IPIV(I)«1 

I1«I*1 

00  11  d»I1,N 
Y-A(d.l) 

IF<Y.LT.O.)Y— Y 

IF(Y.LE.X)G0  TO  11 
X-Y 

IPIV(I) 

11  CONTINUE 

IF(IPIvn). €0.1)00  TO  14 
K*1PI V( I ) 

00  12  d-I.N 
X«A( I , d) 

A(I.d)*A(K.d) 

12  A(K,d)*X 

14  DO  10  d*l  1  ,N 

X— A(d,I)/A(l.I) 

A(d.l)-X 
DO  10  K»I1.N 

10  A(d.X)*A(d,KH'A(I.K) 
RETURN 

I  NO 


1 


SUBROUTINE  BACK(ND.N.IPIV.A.V) 
IMPLICIT  REAL*8( A-H.O-Z ) 
DIMENSION  IPIV(ND) , A(NO.NO) ,V(NO) 
N1»N-1 
DO  10  1*1. N1 
11*1*1 
K*1PIV( I ) 

IF (K . EO . I )G0  TO  11 
X*V<1) 

V( I )*V(K) 

V(K)«X 

11  DO  10  J-11.N 
10  V(d)*V( J)*A( J. 1 )*V( I  ) 

V(N)*V(N)/A(N.N) 

•  DO  IS  II-2.N 
I-N* 1-1 I 
11*1*1 

DO  1C  0*1 1.N 

1C  V<I)*V(I)-A(I.d)*V(J) 

15  V( I )*V( I )/A( I . I ) 

RETURN 

END 

C 

C 


Subroutine  FUN  (Y,F):  User  supplied  subprogram  for  defining  the  dif¬ 
ferential  equations.  F  is  the  vector  of  the  right 
hand  side  of  the  differential  equations.  (See 
examples.) 

Subroutine  DFUN  (Y ,  DF) :  User  supplied  subroutine  for  evaluation  of  the 

Jacobian  of  the  differential  equations.  DF  is  the 
Jacobian  matrix  with  elements  DF  (I,J)*F(i)/Y(J) . 
(See  examples.) 

Subroutine  OUT  (X,  Y,  IH,  Q) :  User  supplied  subprogram  for  output. 

X:  Current  value  of  time. 

Y:  Current  value  of  concentration  vector  at  each 
collocation  point. 

IH:  Number  of  bisections  (unsuccessful  integrations) 
in  the  current  integration  step. 

Q:  Steplength  acceleration  factor. 

STIFF3  also  uses  three  internal  subroutines:  SIRK3,  LU,  and  BACK.  No 
extra  programming  allowances  need  be  made  for  their  use. 

The  first  set  of  listings  following  are  MAIN,  FUN,  OUT,  and  DFUN  for 
the  catalytic  mechanism  including  npectrochemical  absorbances  as  de¬ 
scribed  in  Section  £,  example  Se.  The  second  set  of  listings  are  for 
spherical  symmetry  simple  electron  transfer  including  the  current  cal¬ 
culation  as  described  in  Section  E,  example  4a. 


1 

2 

3 

4 
8 
6 
7 
B 
• 

10 

11 

13 

13 

14 

15 

16 
17 
IS 

19 

20 
21 
22 
23 
;4 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 
46 
49 
80 
81 
83 
83 
64 

85 

86 
87 
68 
69 
•0 


C 

c 

c 

c 

C  ••  PROGRAM  TO  SIMULATE  THE  STANDARD  CATALYTIC  MECHANISM 

C  ••  THE  CONCENTRATION  PROFILES.  CURRENT.  AND 

C  ••  SPECTRCELECTR0CHEM1CAL  absorbances  are  given  ON  OUTPUT. 

c 

c 

c 

IMPLICIT  REAL»6(A-H,0-2) 

DIMENSION  D1F1(30).D1F2(30).DIF3(20).VECT(30).F(30).FOLD(30) 
DIMENSION  V0LD1 ( 30) , YA<  30  > . VK 1 (30) . YK2( 30) . YK3( 30) 

DIMENSION  DF(30.30),DFOLD(30.30).W(30).V(30).VOLD(30) 

DIMENSION  VV(30).YYY(30) 

COMMON  AAA(30.30).8BB(30.30).ROOT(30).WK(30)  .MANUAL 
EXTERNAL  FUN.OUT ,DFUN 
C  ••  SET  NUMBER  OF  SPLINES 
00  1000  MANUAL >1.10 
C  —  GIVE  ORDER  OF  POLYNOMIAL  DESIRED 
N-6 

C  STIPULATE  WHETHER  0  AND  1  ARE  TO  BE  USED  AS  INTERPOLATION  POINTS 
NO*  1 
N1-1 

C  ••  ENTER  ALPHA  ANO  BETA  DETERMINING  POLYNOMIAL  TYPE 
C  ••  (WE  USE  LEGENDRE  HERE) 

AL-O. 

BE-O. 

C  ••  ENTER  PROGRAM  DIMENSION 

N0*30 

C  ••  CALCULATE  ROOTS  OF  SPECIFIED  POLYNOMIAL(COLLOCAT10N  POINTS) 

CALL  dCOBl (ND.N.N0.N1 ,AL .BE ,01F1 .OIF2.OIF3.ROOT ) 
WRITE(6,2)(R00T(1 ). 1*1,6) 

2  F0RMAT(4F 15.12) 

C  ••  CALCULATE  A  MATRIX  OESCRETIZATION  COEFFICIENTS 

DO  4  1-1.8 

CALL  DFOPR(NO. N.N0.N1. 1.1, DIF 1.DIF2.DIF3. ROOT. VECT) 

DO  4  d-1,8 

4  AAA( 1 , d)-VECT(d) 

C  ••  CALCULATE  8  MATRIX  OESCRETIZATION  COEFFICIENTS 
DO  5  1-1,8 

CALL  OFOPR(NO.N,NO.N1.I.2.D!Ft,OIF2.DIF3,RODT.V€CT) 

DO  5  d>1.8 

8  BBB( l.d)-VECT(d) 

C  --  CALCULATE  QUADRATURE  INTEGRATION  WEIGHT  COEFFICIENTS 
DO  791  1-1,8 

CALL  0F0PR(N0.N.N0.N1,I.3.01F1.DIF2.DIF3. ROOT. VECT) 
WKd)-VECT(I) 

791  CONTINUE 

C  ••  ENTER  TIME  LIMITS  BETWEEN  WHICH  INTEGRATION  IS  TO  BE  PERFORMED 
XO-O.DO 
XI -20. DO 

C  ••  ENTER  INTEGRATION  ERROR  TOLERANCE 
BPS- 1.0-06 

C  ••  ENTER  INITIAL  CONDITIONS  FOR  ALL  SPECIES 
DO  300  1-1,6 

vm-t.oo 

V(1«6)-0.D0 

C  ••  SET  WEIGHT  FACTORS  FOR  ANSWERS 
V(!«6)-1.D0 


81  300  W(I)-1.00 

B3  C  —  SET  INITIAL  VALUE  FOR  INTEGRATION  INTERVAL  STEP 

63  H0-1.0-06 

64  C  ••  BEGIN  INTEGRATION 

«*  CALL  STIFF3( 12. 30. 10. FUN. OFUN, OUT, X0.X1.H0.EPS. W.Y, YOLO.  Y0L01.IP. 

««  1 V A , YX 1 , YX2 . VK3 . OF . OFOLO . T , FOLD ) 

67  1000  CONTINUE 

66  STOP 

69  (NO 


1 

SUBROUTINE  FUN(Y.F) 

3 

IMPLICIT  REAL»8(*-H.0-Z> 

3 

DIMENSION  V( 16) ,F( 16) 

4 

COMMON  4*00.30). 88(30. 30). 800(30). WK(30>  .1 

S 

TR-MANUAL 

6 

C  •• 

SET  UP  SPLINE  POINT  MUL' 1PLIER 

7 

VAR* ( 1 . D0/( ( . 100*TR)**2 ) ) 

• 

C  •• 

STATE  THE  0ESCRETI2ED  DIFFERENTIAL  EQUATIONS 

• 

00  350  1-1,6 

10 

T0DER1«0. 

11 

TD0ER2-O. 

13 

00  351  U* 1 .6 

13 

TDDER1*T0DER1*VAR*(BB(I*1.U*1)*V(U>) 

14 

T00ER2-TDDER2*VAR* ( BB ( I* 1 . J* 1 )« V( U*6 ) ) 

IS 

351 

CONTINUE 

16 

F(I)-T0DER1«VAR*BB(I*1.8)*Y(I*6) 

17 

F( I*6)*TDDER2*VAR*BB( 1*1,1 )“Y( 1*6) 

18 

350 

CONTINUE 

18 

RETURN 

30 

END 

31 

C 

33 

c 

1  SUBROUTINE  OUT ( X , V , IH.O) 

2  IMPLICIT  RE4L*B(4-H,0-Z ) 

3  DIMENSION  YY<30)  .YYY(30) 

4  COMMON  44(30,30)  ,B8( 30,30) ,R00( 30) ,WK(30)  .M4NU4L 

5  DIMENSION  V(3C) 

C  WR1TE(6, 130)X,(Y(I),l*1,6),(Y(I+6),I«1.6) 

7  13C  F0RM4TC  TIME  •' ,E20. 15 ,/,4X . '0.00000000000  3(F 16 . 1 1 . IX ) ./ . IX ,3 

0  1(F16. 11. IX ),3X.' 1.00000000000  \//.4X,' 1.00000000000’, 3(F 16. 11. IX). 

8  2/. IX. 3<F16. 11. IX). 3X. '0.00000000000' ) 

10  C  ••  C4LCUL4TE  THE  46S0RB4NCES  FROM  THE  NOW  KNOWN  C0NCENTR4T10N 

11  C  ••  4T  THE  C0LL0C4TI0N  POINTS  (V),  4ND  THE  INTEGR4TI0N 

12  C  ••  COEFFICIENTS  <WK). 

13  G4USS4-0.D0 

14  G4USSB-0.D0 

16  VV( 1 )*0.00 

16  YY(S)*1.D0 

17  YVY ( 1 )• 1 .00 

IB  YYY(8)-0.D0 

18  DO  782  1*1,6 

20  YY(I*1)«Y(I) 

21  YYY(I41)»Y(I4«) 

22  782  CONTINUE 

23  00  783  1*1.8 

24  G4USS4*G4USS44WK(1)*VY(I) 

25  G4US5B*G4U5SB4WK( I )*YYY(I ) 

26  783  CONTINUE 

27  G4U$S4*G4USS4* 1 .0-6* . 100»M4NU4L* 1 .0-3** .5 

2B  G4USSB«G4USSB* 1 .0-6* . 100*M4NU4L* 1 .D-3** .5 

29  WRITE(6.784)G4USS4,G4USSB 

30  784  F0RM4TC  4BS0RB4NCE  4  *'.E20.15.' 

31  RETURN 

32  END 


4BS0RB4NCE  B  •' ,E20. 15.//) 


IMPLICIT  RrAl  *!)<A-  ll.n-r> 
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